Changeset 2461 for trunk/LMDZ.MARS/libf/aeronomars
- Timestamp:
- Feb 16, 2021, 12:31:33 PM (4 years ago)
- Location:
- trunk/LMDZ.MARS/libf/aeronomars
- Files:
-
- 6 edited
-
calchim_mod.F90 (modified) (10 diffs)
-
inichim_newstart.F90 (modified) (2 diffs)
-
moldiff_red.F90 (modified) (10 diffs)
-
photochemistry.F90 (modified) (38 diffs)
-
photolysis_mod.F90 (modified) (1 diff)
-
photolysis_online.F (modified) (5 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/LMDZ.MARS/libf/aeronomars/calchim_mod.F90
r2433 r2461 23 23 igcm_noplus, igcm_n2plus, igcm_hplus, & 24 24 igcm_hco2plus, igcm_hcoplus, igcm_h2oplus, & 25 igcm_h3oplus, igcm_ohplus, igcm_elec, mmol 25 igcm_h3oplus, igcm_ohplus, igcm_elec, & 26 igcm_hdo_vap, igcm_od, igcm_d, igcm_hd, & 27 igcm_do2, igcm_hdo2, mmol 26 28 27 29 use conc_mod, only: mmean ! mean molecular mass of the atmosphere … … 147 149 integer,save :: i_ohplus=0 148 150 integer,save :: i_elec=0 151 integer,save :: i_hdo=0 152 integer,save :: i_od=0 153 integer,save :: i_d=0 154 integer,save :: i_hd=0 155 integer,save :: i_do2=0 156 integer,save :: i_hdo2=0 149 157 150 158 integer :: ig_vl1 … … 154 162 integer :: nquench ! number of quenching + heterogeneous reactions 155 163 integer :: nphotion ! number of photoionizations 164 integer :: nb_reaction_4_ion ! quadratic reactions for ionosphere 165 integer :: nb_reaction_4_deut ! quadratic reactions for deuterium chem 156 166 integer :: nb_phot_max ! total number of photolysis+photoionizations+quenching reactions 157 167 … … 169 179 logical, save :: depos ! switch for dry deposition 170 180 logical, save :: ionchem ! switch for ionospheric chemistry 181 logical, save :: deutchem ! switch for deuterium chemistry 171 182 logical, save :: jonline ! switch for online photodissociation rates or lookup table 172 183 logical, save :: unichim ! only one unified chemical scheme at all … … 217 228 jonline = .true. ! true : on-line calculation of photodissociation rates ! false : lookup table 218 229 ionchem = .false. ! switch for ionospheric chemistry 230 deutchem= .false. ! switch for deuterium chemistry 219 231 depos = .false. ! switch for dry deposition 220 output = . false. ! true : triggers writing of specific outputs (photolysis rates, emission rates, etc)221 222 if (photochem) then223 if (jonline) then224 print*,'calchim: Read UV absorption cross-sections'225 call init_photolysis ! on-line photolysis226 else227 print*,'calchim: Read photolysis lookup table'228 call read_phototable ! off-line photolysis229 end if230 end if231 232 if(.not.unichim) then232 output = .true. ! true : triggers writing of specific outputs (photolysis rates, emission rates, etc) 233 234 ! if (photochem) then 235 ! if (jonline) then 236 ! print*,'calchim: Read UV absorption cross-sections' 237 ! call init_photolysis ! on-line photolysis 238 ! else 239 ! print*,'calchim: Read photolysis lookup table' 240 ! call read_phototable ! off-line photolysis 241 ! end if 242 ! end if 243 244 ! if(.not.unichim) then 233 245 !Read reaction rates from external file if the upper atmospheric 234 246 !chemistry is called 235 call chemthermos_readini236 endif247 ! call chemthermos_readini 248 ! endif 237 249 238 250 ! find index of chemical tracers to use … … 393 405 niq(nbq) = i_h2o 394 406 end if 407 i_hdo=igcm_hdo_vap 408 if (i_hdo == 0) then 409 write(*,*) "calchim: no HDO tracer !!!" 410 ! call abort_physic("calchim","missing hdo_vap tracer",1) 411 write(*,*) "No deuterium chemistry considered" 412 else 413 nbq = nbq + 1 414 niq(nbq) = i_hdo 415 deutchem = .true. 416 write(*,*) "calchim: HDO tracer found in traceur.def" 417 write(*,*) "Deuterium chemistry included" 418 end if 419 i_od=igcm_od 420 if(deutchem) then 421 if (i_od == 0) then 422 write(*,*) "calchim: Error, no OD tracer !!!" 423 write(*,*) "OD is needed if HDO is in traceur.def" 424 call abort_physic("calchim","missing od tracer",1) 425 else 426 nbq = nbq + 1 427 niq(nbq) = i_od 428 end if 429 else 430 if (i_oplus /= 0) then 431 write(*,*) "calchim: Error: OD is present, but HDO is not!!!" 432 write(*,*) "Both must be in traceur.def if deuterium chemistry wanted" 433 call abort_physic("calchim","missing hdo tracer",1) 434 endif 435 endif 436 i_d=igcm_d 437 if(deutchem) then 438 if (i_d == 0) then 439 write(*,*) "calchim: Error, no D tracer !!!" 440 write(*,*) "D is needed if HDO is in traceur.def" 441 call abort_physic("calchim","missing d tracer",1) 442 else 443 nbq = nbq + 1 444 niq(nbq) = i_d 445 end if 446 else 447 if (i_d /= 0) then 448 write(*,*) "calchim: Error: D is present, but HDO is not!!!" 449 write(*,*) "Both must be in traceur.def if deuterium chemistry wanted" 450 call abort_physic("calchim","missing hdo tracer",1) 451 endif 452 endif 453 i_hd=igcm_hd 454 if(deutchem) then 455 if (i_hd == 0) then 456 write(*,*) "calchim: Error, no HD tracer !!!" 457 write(*,*) "HD is needed if HDO is in traceur.def" 458 call abort_physic("calchim","missing hd tracer",1) 459 else 460 nbq = nbq + 1 461 niq(nbq) = i_hd 462 end if 463 else 464 if (i_hd /= 0) then 465 write(*,*) "calchim: Error: HD is present, but HDO is not!!!" 466 write(*,*) "Both must be in traceur.def if deuterium chemistry wanted" 467 call abort_physic("calchim","missing hd tracer",1) 468 endif 469 endif 470 i_do2=igcm_do2 471 if(deutchem) then 472 if (i_do2 == 0) then 473 write(*,*) "calchim: Error, no DO2 tracer !!!" 474 write(*,*) "DO2 is needed if HDO is in traceur.def" 475 call abort_physic("calchim","missing do2 tracer",1) 476 else 477 nbq = nbq + 1 478 niq(nbq) = i_do2 479 end if 480 else 481 if (i_do2 /= 0) then 482 write(*,*) "calchim: Error: DO2 is present, but HDO is not!!!" 483 write(*,*) "Both must be in traceur.def if deuterium chemistry wanted" 484 call abort_physic("calchim","missing do2 tracer",1) 485 endif 486 endif 487 i_hdo2=igcm_hdo2 488 if(deutchem) then 489 if (i_hdo2 == 0) then 490 write(*,*) "calchim: Error, no HDO2 tracer !!!" 491 write(*,*) "HDO2 is needed if HDO is in traceur.def" 492 call abort_physic("calchim","missing hdo2 tracer",1) 493 else 494 nbq = nbq + 1 495 niq(nbq) = i_hdo2 496 end if 497 else 498 if (i_hdo2 /= 0) then 499 write(*,*) "calchim: Error: HDO2 is present, but HDO is not!!!" 500 write(*,*) "Both must be in traceur.def if deuterium chemistry wanted" 501 call abort_physic("calchim","missing hdo2 tracer",1) 502 endif 503 endif 395 504 i_o2plus = igcm_o2plus 396 505 if (i_o2plus == 0) then … … 631 740 endif 632 741 write(*,*) 'calchim: found nbq = ',nbq,' tracers' 633 742 743 write(*,*) 'calchim: tracer indices=',niq(:) 744 745 746 if (photochem) then 747 if (jonline) then 748 print*,'calchim: Read UV absorption cross-sections' 749 !Add two photodissociations in deuterium chemistry included 750 if(deutchem) nphot = nphot + 2 751 call init_photolysis ! on-line photolysis 752 else 753 print*,'calchim: Read photolysis lookup table' 754 call read_phototable ! off-line photolysis 755 end if 756 end if 757 758 if(.not.unichim) then 759 !Read reaction rates from external file if the upper atmospheric 760 !chemistry is called 761 call chemthermos_readini 762 endif 763 764 ! stop 634 765 firstcall = .false. 635 766 end if ! if (firstcall) … … 682 813 lswitch = nlayer + 1 683 814 else 684 if (foundswitch == 0 .and. pplay(ig,l) < 1 .e-1) then815 if (foundswitch == 0 .and. pplay(ig,l) < 10.) then 685 816 lswitch = l 686 817 foundswitch = 1 … … 698 829 if (photochem) then 699 830 ! set number of reactions, depending on ion chemistry or not 831 nb_reaction_4_ion = 64 832 nb_reaction_4_deut = 35 833 834 !Default numbers if no ion and no deuterium chemistry included 835 836 nb_reaction_4_max = 31 ! set number of bimolecular reactions 837 nb_reaction_3_max = 6 ! set number of quadratic reactions 838 nquench = 9 ! set number of quenching + heterogeneous 839 nphotion = 0 ! set number of photoionizations 840 700 841 if (ionchem) then 701 nb_reaction_4_max = 95 ! set number of bimolecular reactions 702 nb_reaction_3_max = 6 ! set number of quadratic reactions 703 nquench = 9 ! set number of quenching + heterogeneous reactions 842 nb_reaction_4_max = nb_reaction_4_max + nb_reaction_4_ion 704 843 nphotion = 18 ! set number of photoionizations 705 else 706 nb_reaction_4_max = 31 ! set number of bimolecular reactions 707 nb_reaction_3_max = 6 ! set number of quadratic reactions 708 nquench = 9 ! set number of quenching + heterogeneous reactions 709 nphotion = 0 ! set number of photoionizations 844 endif 845 if(deutchem) then 846 nb_reaction_4_max = nb_reaction_4_max + nb_reaction_4_deut 710 847 end if 711 848 … … 717 854 ! call main photochemistry routine 718 855 719 call photochemistry(nlayer,nq,nbq,ionchem,nb_reaction_3_max, & 720 nb_reaction_4_max,nphot,nb_phot_max,nphotion, & 856 call photochemistry(nlayer,nq,nbq,ionchem,deutchem, & 857 nb_reaction_3_max,nb_reaction_4_max,nphot, & 858 nb_phot_max,nphotion, & 721 859 jonline,ig,lswitch,zycol,szacol,ptimestep, & 722 860 zpress,zlocal,ztemp,ztelec,zdens,zmmean, & -
trunk/LMDZ.MARS/libf/aeronomars/inichim_newstart.F90
r2409 r2461 108 108 igcm_n2d=0 109 109 igcm_he=0 110 igcm_hdo_vap=0 111 igcm_od=0 112 igcm_d=0 113 igcm_hd=0 114 igcm_do2=0 115 igcm_hdo2=0 110 116 igcm_co2plus=0 111 117 igcm_oplus=0 … … 277 283 igcm_h2o_ice = iq 278 284 mmol(igcm_h2o_ice) = 18. 285 count = count + 1 286 end if 287 if (noms(iq) == "hdo_vap") then 288 igcm_hdo_vap = iq 289 mmol(igcm_hdo_vap) = 19. 290 count = count + 1 291 end if 292 if (noms(iq) == "od") then 293 igcm_od = iq 294 mmol(igcm_od) = 18. 295 count = count + 1 296 end if 297 if (noms(iq) == "d") then 298 igcm_d = iq 299 mmol(igcm_d) = 2. 300 count = count + 1 301 end if 302 if (noms(iq) == "hd") then 303 igcm_d = iq 304 mmol(igcm_d) = 3. 305 count = count + 1 306 end if 307 if (noms(iq) == "do2") then 308 igcm_do2 = iq 309 mmol(igcm_do2) = 34. 310 count = count + 1 311 end if 312 if (noms(iq) == "hdo2") then 313 igcm_hdo2 = iq 314 mmol(igcm_hdo2) = 35. 279 315 count = count + 1 280 316 end if -
trunk/LMDZ.MARS/libf/aeronomars/moldiff_red.F90
r2157 r2461 43 43 real*8 masseU,kBolt,RgazP,Rmars,Grav,Mmars 44 44 real*8 rho0,D0,T0,H0,time0,dZ,time,dZraf,tdiff,Zmin,Zmax 45 real*8 FacEsc,invsgmu,PhiEscH,PhiEscH2 45 real*8 FacEsc,invsgmu,PhiEscH,PhiEscH2,PhiEscD 46 46 real*8 hp(nlayer) 47 47 real*8 pp(nlayer) … … 64 64 real*8,dimension(:),allocatable,save :: wi,Wad,Uthermal,Lambdaexo,Hspecie 65 65 real*8,dimension(:),allocatable,save :: Mtot1,Mtot2,Mraf1,Mraf2 66 integer,parameter :: ListeDiffNb=1 566 integer,parameter :: ListeDiffNb=16 67 67 character(len=20),dimension(ListeDiffNb) :: ListeDiff 68 68 real*8,parameter :: pi=3.141592653 … … 75 75 integer,save :: i_h2 76 76 integer,save :: i_h 77 integer,save :: i_d 77 78 ! vertical index limit for the molecular diffusion 78 79 integer,save :: il0 … … 130 131 ListeDiff(14)='n' 131 132 ListeDiff(15)='he' 133 ListeDiff(16)='hdo_vap' 132 134 i_h=1000 133 135 i_h2=1000 136 i_d=1000 134 137 ! On regarde quelles especes diffusables sont presentes 135 138 … … 162 165 if (noms(nn) .eq. 'h') i_h=n 163 166 if (noms(nn) .eq. 'h2') i_h2=n 167 if (noms(nn) .eq. 'd') i_d=n 164 168 endif 165 169 enddo … … 222 226 PhiEscH=0D0 223 227 PhiEscH2=0D0 228 PhiEscD=0D0 224 229 225 230 do ig=1,ngrid … … 403 408 & (Rmars+Zraf(nlraf))/kbolt/Traf(nlraf) 404 409 wi(nn)=0D0 405 if (nn .eq. i_h .or. nn .eq. i_h2 ) then410 if (nn .eq. i_h .or. nn .eq. i_h2 .or. nn .eq. i_d) then 406 411 wi(nn)=Uthermal(nn)/2d0/sqrt(PI)*exp(-Lambdaexo(nn))* & 407 412 & (Lambdaexo(nn)+1d0) … … 409 414 enddo 410 415 411 ! print*,'wi',wi(i_h),wi(i_h2), Uthermal,Lambdaexo,mmol(gcmind(:))416 ! print*,'wi',wi(i_h),wi(i_h2),wi(i_d),Uthermal,Lambdaexo,mmol(gcmind(:)) 412 417 ! print*,'wi',wi 413 ! stop414 418 415 419 ! Compute time step for diffusion … … 594 598 ! the trend only at the end 595 599 596 PhiEscH=PhiEscH+wi(i_h)*Nrafk(nlraf,i_h)*cell_area(ig) ! in s-1 597 PhiEscH2=PhiEscH2+wi(i_h2)*Nrafk(nlraf,i_h2)*cell_area(ig) ! in s-1 (U in m/s, aire in m2, Nrafk in m-3) 598 ! print*,'test',ig,wi(i_h),Nrafk(nlraf,i_h),wi(i_h2),Nrafk(nlraf,i_h2),cell_area(ig),PhiEscH,PhiEscH2,i_h,i_h2 600 if (i_h .ne. 1000) PhiEscH=PhiEscH+wi(i_h)*Nrafk(nlraf,i_h)*cell_area(ig) ! in s-1 601 if (i_h2 .ne. 1000) PhiEscH2=PhiEscH2+wi(i_h2)*Nrafk(nlraf,i_h2)*cell_area(ig) ! in s-1 (U in m/s, aire in m2, Nrafk in m-3) 602 if (i_d .ne. 1000) PhiEscD=PhiEscD+wi(i_d)*Nrafk(nlraf,i_d)*cell_area(ig) 603 ! print*,'test',ig,wi(i_h),wi(i_h2),Nrafk(nlraf,i_h),Nrafk(nlraf,i_h2),Nrafk(nlraf,i_d),cell_area(ig),PhiEscH,PhiEscH2,i_h,i_h2,i_d,PhiEscD 599 604 ! stop 600 605 … … 651 656 652 657 enddo ! ig loop 653 ! print*,'Escape flux H, H2 (s-1)',PhiEscH,PhiEscH2658 ! print*,'Escape flux H, H2,D (s-1)',PhiEscH,PhiEscH2,PhiEscD 654 659 655 660 return -
trunk/LMDZ.MARS/libf/aeronomars/photochemistry.F90
r2433 r2461 13 13 !***************************************************************** 14 14 15 subroutine photochemistry(nlayer, nq, nesp, ionchem, nb_reaction_3_max, & 16 nb_reaction_4_max, nphot, nb_phot_max, nphotion, & 15 subroutine photochemistry(nlayer, nq, nesp, ionchem, deutchem, & 16 nb_reaction_3_max, nb_reaction_4_max, nphot, & 17 nb_phot_max, nphotion, & 17 18 jonline, ig, lswitch, zycol, sza, ptimestep, press, & 18 19 alt, temp, temp_elect, dens, zmmean, & … … 45 46 ! number of photoionizations 46 47 logical, intent(in) :: ionchem! switch for ion chemistry 48 logical, intent(in) :: deutchem 49 ! switch for deuterium chemistry 47 50 logical, intent(in) :: jonline! switch for on-line calculation of photolysis rates 48 51 integer :: ig ! grid point index … … 89 92 90 93 ! tracer indexes in the chemistry: 91 94 ! Species always considered in the chemistry 92 95 integer,parameter :: i_co2 = 1 93 96 integer,parameter :: i_co = 2 … … 107 110 integer,parameter :: i_no2 = 16 108 111 integer,parameter :: i_n2 = 17 109 integer,parameter :: i_co2plus = 18 110 integer,parameter :: i_oplus = 19 111 integer,parameter :: i_o2plus = 20 112 integer,parameter :: i_noplus = 21 113 integer,parameter :: i_coplus = 22 114 integer,parameter :: i_cplus = 23 115 integer,parameter :: i_n2plus = 24 116 integer,parameter :: i_nplus = 25 117 integer,parameter :: i_hplus = 26 118 integer,parameter :: i_hco2plus= 27 119 integer,parameter :: i_hcoplus = 28 120 integer,parameter :: i_h2oplus = 29 121 integer,parameter :: i_h3oplus = 30 122 integer,parameter :: i_ohplus = 31 123 integer,parameter :: i_elec = 32 112 !Species that may or not be included 113 !Their indices will be defined, if species are to be considered, in firstcall 114 integer,save :: i_co2plus = 0 115 integer,save :: i_oplus = 0 116 integer,save :: i_o2plus = 0 117 integer,save :: i_noplus = 0 118 integer,save :: i_coplus = 0 119 integer,save :: i_cplus = 0 120 integer,save :: i_n2plus = 0 121 integer,save :: i_nplus = 0 122 integer,save :: i_hplus = 0 123 integer,save :: i_hco2plus = 0 124 integer,save :: i_hcoplus = 0 125 integer,save :: i_h2oplus = 0 126 integer,save :: i_h3oplus = 0 127 integer,save :: i_ohplus = 0 128 integer,save :: i_elec = 0 129 integer,save :: i_hdo = 0 130 integer,save :: i_od = 0 131 integer,save :: i_d = 0 132 integer,save :: i_hd = 0 133 integer,save :: i_do2 = 0 134 integer,save :: i_hdo2 = 0 135 136 !integer,parameter :: i_co2plus = 18 137 !integer,parameter :: i_oplus = 19 138 !integer,parameter :: i_o2plus = 20 139 !integer,parameter :: i_noplus = 21 140 !integer,parameter :: i_coplus = 22 141 !integer,parameter :: i_cplus = 23 142 !integer,parameter :: i_n2plus = 24 143 !integer,parameter :: i_nplus = 25 144 !integer,parameter :: i_hplus = 26 145 !integer,parameter :: i_hco2plus= 27 146 !integer,parameter :: i_hcoplus = 28 147 !integer,parameter :: i_h2oplus = 29 148 !integer,parameter :: i_h3oplus = 30 149 !integer,parameter :: i_ohplus = 31 150 !integer,parameter :: i_elec = 32 151 !integer,parameter :: i_hdo = 33 152 !integer,parameter :: i_od = 34 153 !integer,parameter :: i_d = 35 154 !integer,parameter :: i_hd = 36 155 !integer,parameter :: i_do2 = 37 156 !integer,parameter :: i_hdo2 = 38 157 158 integer :: i_last 124 159 125 160 !Tracer indexes for photionization coeffs … … 173 208 174 209 if (firstcall) then 210 !Indices for species that may or not be considered 211 i_last = i_n2 212 if(ionchem) then 213 i_co2plus = i_last + 1 214 i_oplus = i_last + 2 215 i_o2plus = i_last + 3 216 i_noplus = i_last + 4 217 i_coplus = i_last + 5 218 i_cplus = i_last + 6 219 i_n2plus = i_last + 7 220 i_nplus = i_last + 8 221 i_hplus = i_last + 9 222 i_hco2plus = i_last + 10 223 i_hcoplus = i_last + 11 224 i_h2oplus = i_last + 12 225 i_h3oplus = i_last + 13 226 i_ohplus = i_last + 14 227 i_elec = i_last + 15 228 !Update last used index 229 i_last = i_elec 230 endif 231 if(deutchem) then 232 i_hdo = i_last + 1 233 i_od = i_last + 2 234 i_d = i_last + 3 235 i_hd = i_last + 4 236 i_do2 = i_last + 5 237 i_hdo2 = i_last + 6 238 i_last = i_hdo2 239 endif 240 print*,'photochemistry: check number of indices',i_last,nesp 175 241 print*,'photochemistry: initialize indexes' 176 242 call indice(nb_reaction_3_max,nb_reaction_4_max, & 177 nb_phot_max, ionchem, i_co2, i_co, i_o, i_o1d, i_o2,&178 i_o 3, i_h,i_h2, i_oh, i_ho2, i_h2o2, i_h2o,&179 i_ n, i_n2d, i_no, i_no2, i_n2, i_co2plus,&243 nb_phot_max, ionchem, deutchem, i_co2, i_co, i_o, & 244 i_o1d, i_o2, i_o3, i_h,i_h2, i_oh, i_ho2, i_h2o2, & 245 i_h2o, i_n, i_n2d, i_no, i_no2, i_n2, i_co2plus, & 180 246 i_oplus, i_o2plus, i_noplus, i_coplus, i_cplus, & 181 247 i_n2plus, i_nplus, i_hplus, i_hco2plus, i_hcoplus, & 182 i_h2oplus, i_h3oplus, i_ohplus, i_elec) 248 i_h2oplus, i_h3oplus, i_ohplus, i_elec, & 249 i_hdo, i_od, i_d, i_hd, i_do2, i_hdo2) 183 250 firstcall = .false. 184 251 end if … … 188 255 !=================================================================== 189 256 190 call gcmtochim(nlayer, ionchem, nq, zycol, lswitch, nesp,&191 i_co2, i_co, i_o, i_o1d, i_o2, i_o3, i_h,&257 call gcmtochim(nlayer, ionchem, deutchem, nq, zycol, lswitch, & 258 nesp, i_co2, i_co, i_o, i_o1d, i_o2, i_o3, i_h, & 192 259 i_h2, i_oh, i_ho2, i_h2o2, i_h2o, & 193 260 i_n, i_n2d, i_no, i_no2, i_n2, i_co2plus, & … … 195 262 i_n2plus, i_nplus, i_hplus, i_hco2plus, & 196 263 i_hcoplus, i_h2oplus, i_h3oplus, i_ohplus, & 197 i_elec, dens, rm, c) 264 i_elec, i_hdo, i_od, i_d, i_hd, i_do2, i_hdo2, & 265 dens, rm, c) 198 266 199 267 !=================================================================== … … 205 273 if (jonline) then 206 274 if (sza <= 113.) then ! day at 300 km 207 call photolysis_online(nlayer, nb_phot_max, alt, press, temp, zmmean, & 275 call photolysis_online(nlayer, deutchem, nb_phot_max, alt, & 276 press, temp, zmmean, & 208 277 i_co2, i_co, i_o, i_o1d, i_o2, i_o3, i_h, & 209 278 i_h2, i_oh, i_ho2, i_h2o2, i_h2o, & … … 306 375 hetero_ice = .false. 307 376 308 call reactionrates(nlayer, ionchem, nb_reaction_3_max, nb_reaction_4_max, & 377 call reactionrates(nlayer, ionchem, deutchem, & 378 nb_reaction_3_max, nb_reaction_4_max, & 309 379 nb_phot_max, nphotion, lswitch, dens, c(:,i_co2), & 310 380 c(:,i_o2), c(:,i_o), c(:,i_n2), press, temp, & … … 429 499 !=================================================================== 430 500 431 call chimtogcm(nlayer, ionchem, nq, zycol, lswitch, nesp,&432 i_co2, i_co, i_o, i_o1d, i_o2, i_o3, i_h,&501 call chimtogcm(nlayer, ionchem, deutchem, nq, zycol, lswitch, & 502 nesp, i_co2, i_co, i_o, i_o1d, i_o2, i_o3, i_h, & 433 503 i_h2, i_oh, i_ho2, i_h2o2, i_h2o, & 434 504 i_n, i_n2d, i_no, i_no2, i_n2, i_co2plus, & … … 436 506 i_n2plus, i_nplus, i_hplus, i_hco2plus, & 437 507 i_hcoplus, i_h2oplus, i_h3oplus, i_ohplus, & 438 i_elec, dens, c) 508 i_elec, i_hdo, i_od, i_d, i_hd, i_do2, i_hdo2, & 509 dens, c) 439 510 contains 440 511 … … 560 631 !====================================================================== 561 632 562 subroutine reactionrates(nlayer, ionchem, nb_reaction_3_max, nb_reaction_4_max, & 633 subroutine reactionrates(nlayer, ionchem, deutchem , & 634 nb_reaction_3_max, nb_reaction_4_max, & 563 635 nb_phot_max, nphotion, lswitch, dens, co2, o2, o, & 564 636 n2, press, t, t_elect, hetero_dust, hetero_ice, & … … 591 663 integer, intent(in) :: nphotion ! number of photoionizations 592 664 logical, intent(in) :: ionchem 665 logical, intent(in) :: deutchem 593 666 integer :: lswitch ! interface level between lower 594 667 ! atmosphere and thermosphere chemistries … … 624 697 real, dimension(nlayer) :: a001, a002, a003, & 625 698 b001, b002, b003, b004, b005, b006, b007, & 626 b008, b009, 699 b008, b009, b010, b011, b012, & 627 700 c001, c002, c003, c004, c005, c006, c007, & 628 701 c008, c009, c010, c011, c012, c013, c014, & … … 631 704 d008, d009, d010, d011, d012, & 632 705 e001, e002, & 706 f001, f002, f003, f004, f005, f006, f007, & 707 f008, f009, f010, f011, f012, f013, f014, & 708 f015, f016, f017, f018, f019, f020, f021, & 709 f022, f023, f024, f025, f026, f027, f028, & 710 f029, f030, f031, f032, & 633 711 i001, i002, i003, i004, i005, i006, & 634 712 i007, i008, i009, i010, i011, i012, & … … 766 844 b009(:) = 1.5e-10*0.05 767 845 846 if(deutchem) then 847 ! 848 !--- b010: o(1d) + hdo -> od + oh 849 b010(:) = b002(:) 850 851 nb_reaction_4 = nb_reaction_4 + 1 852 v_4(:,nb_reaction_4) = b010(:) 853 854 ! 855 !--- b011: o(1d) + hd -> h + od 856 857 !Laurent et al., 1995 858 859 b011(:) = 1.3e-10 860 861 nb_reaction_4 = nb_reaction_4 + 1 862 v_4(:,nb_reaction_4) = b011(:) 863 864 ! 865 !--- b012: o(1d) + hd -> d + oh 866 867 !Laurent et al., 1995 868 869 b012(:) = 1.0e-10 870 871 nb_reaction_4 = nb_reaction_4 + 1 872 v_4(:,nb_reaction_4) = b012(:) 873 874 endif !deutchem 768 875 !---------------------------------------------------------------------- 769 876 ! hydrogen reactions … … 963 1070 nb_reaction_3 = nb_reaction_3 + 1 964 1071 v_3(:,nb_reaction_3) = c018(:) 1072 965 1073 966 1074 !---------------------------------------------------------------------- … … 1137 1245 nb_reaction_4 = nb_reaction_4 + 1 1138 1246 v_4(:,nb_reaction_4) = e002(:) 1247 1248 1249 !---------------------------------------------------------------------- 1250 ! deuterium reactions 1251 ! only if deutchem = true 1252 !---------------------------------------------------------------------- 1253 1254 if(deutchem) then 1255 1256 !--- f001: od + oh -> hdo + o 1257 1258 ! Bedjanian, Y.; Le Bras, G.; and Poulet, G., 1259 !Kinetic Study of OH + OH and OD + OD Reactions, 1260 !J. Phys. Chem. A, 103, pp. 7017 - 7025, 1999 1261 1262 f001(:) = 1.5e-12 1263 1264 nb_reaction_4 = nb_reaction_4 + 1 1265 v_4(:,nb_reaction_4) = f001(:) 1266 1267 !--- f002: od + h2 -> HDO + H 1268 1269 !Talukdar, R.K. et al., 1270 !Kinetics of hydroxyl radical reactions with isotopically labeled hydrogen, 1271 !J. Phys. Chem., 100, pp. 3037 - 3043, 1996 1272 1273 f002(:) = 7.41e-15 1274 1275 nb_reaction_4 = nb_reaction_4 + 1 1276 v_4(:,nb_reaction_4) = f002(:) 1277 1278 !--- f003: od + ho2 -> hdo + o2 1279 1280 f003(:) = c007(:) 1281 1282 nb_reaction_4 = nb_reaction_4 + 1 1283 v_4(:,nb_reaction_4) = f003(:) 1284 1285 !--- f004: od + h2o2 -> hdo + ho2 1286 1287 !Vaghjiani, G.L. & Ravishankara, A.R., 1288 !Reactions of OH and OD with H2O2 and D2O2, 1289 !J. Phys. Chem., 93, pp. 7833 - 7837, 1989; 1290 1291 f004(:) = 1.79e-12 1292 1293 nb_reaction_4 = nb_reaction_4 + 1 1294 v_4(:,nb_reaction_4) = f004(:) 1295 1296 !--- f005: o + od -> o2 + d 1297 1298 !Following Yung+1998, rate equal to that of O + OH -> O2 + H (c002) 1299 1300 f005(:) = c002(:) 1301 1302 nb_reaction_4 = nb_reaction_4 + 1 1303 v_4(:,nb_reaction_4) = f005(:) 1304 1305 !--- f006: od + h2 -> h2o + d 1306 1307 !Following Yung+1998, rate equal to that of OH + H2 -> H2O + H (c010) 1308 1309 f006(:) = c010(:) 1310 1311 nb_reaction_4 = nb_reaction_4 + 1 1312 v_4(:,nb_reaction_4) = f006(:) 1313 1314 !--- f007: od + h -> oh + d 1315 1316 !Rate following Yung+1988 and the rate of the inverse reaction (f012) from Atahan et al. J. Chem. Phys. 2005 1317 1318 f007(:) = 1.61e-10*((298./t(:))**0.32)*exp(-16./t(:)) 1319 1320 nb_reaction_4 = nb_reaction_4 + 1 1321 v_4(:,nb_reaction_4) = f007(:) 1322 1323 !--- f008: co + od -> co2 + d 1324 1325 !Following Yung+1988 rate equal to that of reaction CO + OH -> CO2 + H 1326 1327 f008(:) = e001(:) 1328 1329 nb_reaction_4 = nb_reaction_4 + 1 1330 v_4(:,nb_reaction_4) = f008(:) 1331 1332 !--- f009: o3 + d -> o2 + od 1333 1334 !Rate from NIST, Yu & Varandas, J. Chem. Soc. Faraday Trans., 93, 2651-2656, 1997 1335 1336 f009(:) = 7.41e-11*exp(-379./t(:)) 1337 1338 nb_reaction_4 = nb_reaction_4 + 1 1339 v_4(:,nb_reaction_4) = f009(:) 1340 1341 !--- f010: HO2 + D -> OH + OD 1342 1343 !Following Yung+1988, rate equal to 0.71 times the rate for reaction HO2 + H -> OH + OH (c004) 1344 1345 f010(:) = 0.71*c004(:) 1346 1347 nb_reaction_4 = nb_reaction_4 + 1 1348 v_4(:,nb_reaction_4) = f010(:) 1349 1350 !--- f011: HO2 + D -> HDO + O 1351 1352 !Following Yung+1988, rate equal to 0.71 times the rate for reaction HO2 + H -> H2O + O (c006) 1353 1354 f011(:) = 0.71*c006(:) 1355 1356 nb_reaction_4 = nb_reaction_4 + 1 1357 v_4(:,nb_reaction_4) = f011(:) 1358 1359 !--- f012: OH + D -> OD + H 1360 1361 !Rate from NIST, Atahan et al. J. Chem. Phys. 2005 1362 1363 f012(:) = 1.16e-10*((298./t(:))**0.32)*exp(-16./t(:)) 1364 1365 nb_reaction_4 = nb_reaction_4 + 1 1366 v_4(:,nb_reaction_4) = f012(:) 1367 1368 !--- f013: h + d + co2 -> hd + co2 1369 1370 !According to Yung et al. 1988, rate equal to that of H + H + CO2 1371 !(reaction c018). Source: baulch et al., 2005 1372 1373 f013(:) = c018(:) 1374 1375 nb_reaction_4 = nb_reaction_4 + 1 1376 v_4(:,nb_reaction_4) = f013(:) 1377 1378 !--- f014: D + HO2 -> HD + O2 1379 1380 !According to Yung et al., rate equal to 0.71 times the rate of 1381 ! H + HO2 -> H2 + O2 (reaction c005, source JPL 2019) 1382 1383 f014(:) = 0.71*c005(:) 1384 1385 nb_reaction_4 = nb_reaction_4 + 1 1386 v_4(:,nb_reaction_4) = f014(:) 1387 1388 !--- f015: OH + HD -> HDO + H 1389 1390 !Talukdar et al., Kinetics of hydroxyl radical reactions with 1391 !isotopically labeled hydrogen, J. Phys. Chem. 100, 3037-3043, 1996 1392 1393 f015(:) = 8.5e-13*exp(-2130./t(:)) 1394 1395 nb_reaction_4 = nb_reaction_4 + 1 1396 v_4(:,nb_reaction_4) = f015(:) 1397 1398 !--- f016: OH + HD -> H2O + D 1399 1400 !Talukdar et al., 1996 1401 1402 f016(:) = 4.15e-12*exp(-2130./t(:)) 1403 1404 nb_reaction_4 = nb_reaction_4 + 1 1405 v_4(:,nb_reaction_4) = f016(:) 1406 1407 !--- f017: D + O2 + CO2 -> DO2 + CO2 1408 1409 !Breshears et al., Room-temperature rate constants for the reaction 1410 !D + O2 + M -> DO2 + M, M=Ar, D2, CO2 and F2. J. Chem. Soc. Faraday 1411 !Trans., 87, 2337-2355 (1991) 1412 1413 f017(:) = 1.59e-31 1414 1415 nb_reaction_4 = nb_reaction_4 + 1 1416 v_4(:,nb_reaction_4) = f017(:) 1417 1418 !--- f018: OD + O3 -> DO2 + O2 1419 1420 !According to Yung+1988, rate equal to that of reaccion 1421 !OH + O3 -> HO2 + O2, (reaction c014) 1422 1423 f018(:) = c014(:) 1424 1425 nb_reaction_4 = nb_reaction_4 + 1 1426 v_4(:,nb_reaction_4) = f018(:) 1427 1428 !--- f019: D + HO2 -> DO2 + H 1429 1430 !Yung et al., 1988 1431 1432 f019(:) = 1.0e-10 1433 1434 nb_reaction_4 = nb_reaction_4 + 1 1435 v_4(:,nb_reaction_4) = f019(:) 1436 1437 !--- f020: O + DO2 -> OD + O2 1438 1439 !According to Yung+1988, rate equal to that of O + HO2 -> OH + O2 1440 ! -> reaction c001 1441 1442 f020(:) = c001(:) 1443 1444 nb_reaction_4 = nb_reaction_4 + 1 1445 v_4(:,nb_reaction_4) = f020(:) 1446 1447 !--- f021: H + DO2 -> OH + OD 1448 1449 !According to Yung+1988, rate equal to that of H + HO2 -> OH + OH 1450 ! -> reaction c004 1451 1452 f021(:) = c004(:) 1453 1454 nb_reaction_4 = nb_reaction_4 + 1 1455 v_4(:,nb_reaction_4) = f021(:) 1456 1457 !--- f022: H + DO2 -> HD + O2 1458 1459 !According to Yung+1988, rate equal to that of H + HO2 -> H2 + O2 1460 ! -> reaction c005 1461 1462 f022(:) = c005(:) 1463 1464 nb_reaction_4 = nb_reaction_4 + 1 1465 v_4(:,nb_reaction_4) = f022(:) 1466 1467 !--- f023: H + DO2 -> HDO + O 1468 1469 !According to Yung+1988, rate equal to that of H + HO2 -> H2O + O 1470 ! -> reaction c006 1471 1472 f023(:) = c006(:) 1473 1474 nb_reaction_4 = nb_reaction_4 + 1 1475 v_4(:,nb_reaction_4) = f023(:) 1476 1477 !--- f024: H + DO2 -> HO2 + D 1478 1479 !Yung+1988 1480 1481 f024(:) = 1.85e-10*exp(-890./t(:)) 1482 1483 nb_reaction_4 = nb_reaction_4 + 1 1484 v_4(:,nb_reaction_4) = f024(:) 1485 1486 !--- f025: OH + DO2 -> HDO + O2 1487 1488 !According to Yung+1988, rate equal to that of OH + HO2 -> H2O + O2 1489 ! -> reaction c007 1490 1491 f025(:) = c007(:) 1492 1493 nb_reaction_4 = nb_reaction_4 + 1 1494 v_4(:,nb_reaction_4) = f025(:) 1495 1496 !--- f026: DO2 + O3 -> OD + O2 + O2 1497 1498 !According to Yung+1988, rate equal to that of the reaction 1499 ! HO2 + O3 -> OH + O2 + O2 -> reaction c015 1500 1501 f026(:) = c015(:) 1502 1503 nb_reaction_4 = nb_reaction_4 + 1 1504 v_4(:,nb_reaction_4) = f026(:) 1505 1506 !--- f027: OD + OH + CO2 -> HDO2 + CO2 1507 1508 !According to Yung+1988, rate equal to that of the reaction 1509 ! OH + OH + CO2 -> H2O2 + CO2 (reaction c017) 1510 1511 f027(:) = c017(:) 1512 1513 nb_reaction_4 = nb_reaction_4 + 1 1514 v_4(:,nb_reaction_4) = f027(:) 1515 1516 !--- f028: DO2 + HO2 -> HDO2 + O2 1517 1518 !According to Yung+1988, rate equal to that of HO2 + HO2 -> H2O2 + O2 1519 ! (reaction c008) 1520 1521 f028(:) = c008(:) 1522 1523 nb_reaction_4 = nb_reaction_4 + 1 1524 v_4(:,nb_reaction_4) = f028(:) 1525 1526 !--- f029: O + HDO2 -> OD + HO2 1527 1528 !According to Yung+1988, rate half that of O + H2O2 -> OH + HO2 1529 ! (reaction c012) 1530 1531 f029(:) = 0.5*c012(:) 1532 1533 nb_reaction_4 = nb_reaction_4 + 1 1534 v_4(:,nb_reaction_4) = f029(:) 1535 1536 !--- f030: O + HDO2 -> OH + DO2 1537 1538 !According to Yung+1988, rate half that of O + H2O2 -> OH + HO2 1539 ! (reaction c012) 1540 1541 f030(:) = 0.5*c012(:) 1542 1543 nb_reaction_4 = nb_reaction_4 + 1 1544 v_4(:,nb_reaction_4) = f030(:) 1545 1546 !--- f031: OH + HDO2 -> HDO + HO2 1547 1548 !According to Yung+1988, rate half that of OH + H2O2 -> H2O + HO2 1549 ! (reaction c009) 1550 1551 f031(:) = 0.5*c009(:) 1552 1553 nb_reaction_4 = nb_reaction_4 + 1 1554 v_4(:,nb_reaction_4) = f031(:) 1555 1556 1557 !--- f032: OH + HDO2 -> H2O + DO2 1558 1559 !According to Yung+1988, rate half that of OH + H2O2 -> H2O + HO2 1560 ! (reaction c009) 1561 1562 f032(:) = 0.5*c009(:) 1563 1564 nb_reaction_4 = nb_reaction_4 + 1 1565 v_4(:,nb_reaction_4) = f032(:) 1566 1567 endif !deutchem 1139 1568 1140 1569 !---------------------------------------------------------------------- … … 1891 2320 1892 2321 subroutine indice(nb_reaction_3_max, nb_reaction_4_max, & 1893 nb_phot_max, ionchem, i_co2, i_co, i_o, i_o1d, i_o2,&1894 i_o 3, i_h,i_h2, i_oh, i_ho2, i_h2o2, i_h2o,&1895 i_ n, i_n2d, i_no, i_no2, i_n2, i_co2plus,&2322 nb_phot_max, ionchem, deutchem, i_co2, i_co, i_o, & 2323 i_o1d, i_o2, i_o3, i_h,i_h2, i_oh, i_ho2, i_h2o2, & 2324 i_h2o, i_n, i_n2d, i_no, i_no2, i_n2, i_co2plus, & 1896 2325 i_oplus, i_o2plus, i_noplus, i_coplus, i_cplus, & 1897 2326 i_n2plus, i_nplus, i_hplus, i_hco2plus, i_hcoplus, & 1898 i_h2oplus, i_h3oplus, i_ohplus, i_elec) 2327 i_h2oplus, i_h3oplus, i_ohplus, i_elec, & 2328 i_hdo, i_od, i_d, i_hd, i_do2, i_hdo2) 1899 2329 1900 2330 !================================================================ … … 1921 2351 i_co2plus, i_oplus, i_o2plus, i_noplus, i_coplus, & 1922 2352 i_cplus, i_n2plus, i_nplus, i_hplus, i_hco2plus, & 1923 i_hcoplus, i_h2oplus, i_h3oplus, i_ohplus, i_elec 2353 i_hcoplus, i_h2oplus, i_h3oplus, i_ohplus, i_elec, & 2354 i_hdo, i_od, i_d, i_hd, i_do2, i_hdo2 1924 2355 integer, intent(in) :: nb_reaction_3_max 1925 2356 ! number of quadratic reactions … … 1929 2360 ! number of processes treated numerically as photodissociations 1930 2361 logical, intent(in) :: ionchem 2362 logical, intent(in) :: deutchem 1931 2363 1932 2364 ! local … … 2049 2481 indice_phot(nb_phot) = z3spec(1.0, i_n2, 1.0, i_n2d, 1.0, i_n) 2050 2482 2051 !=========================================================== 2052 ! HDO + hv -> products 2053 !=========================================================== 2054 2055 !nb_phot = nb_phot + 1 2056 2057 !indice_phot(nb_phot) = z3spec(1.0, i_hdo, 0.0, i_dummy, 0.0, i_dummy) 2483 !Only if deuterium chemistry included 2484 if (deutchem) then 2485 !=========================================================== 2486 ! HDO + hv -> OD + H 2487 !=========================================================== 2488 2489 nb_phot = nb_phot + 1 2490 2491 indice_phot(nb_phot) = z3spec(1.0, i_hdo, 1.0, i_h, 1.0, i_od) 2492 2493 !=========================================================== 2494 ! HDO + hv -> D + OH 2495 !=========================================================== 2496 2497 nb_phot = nb_phot + 1 2498 2499 indice_phot(nb_phot) = z3spec(1.0, i_hdo, 1.0, i_d, 1.0, i_oh) 2500 2501 endif !deutchem 2058 2502 2059 2503 !Only if ion chemistry included … … 2273 2717 indice_4(nb_reaction_4) = z4spec(1.0, i_o1d, 1.0, i_o3, 1.0, i_o2, 2.0, i_o) 2274 2718 2719 2720 if(deutchem) then 2721 !=========================================================== 2722 ! b010 : O(1D) + HDO -> OD + OH 2723 !=========================================================== 2724 2725 nb_reaction_4 = nb_reaction_4 + 1 2726 2727 indice_4(nb_reaction_4) = z4spec(1.0, i_o1d, 1.0, i_hdo, 1.0, i_od, 1.0, i_oh) 2728 2729 2730 !=========================================================== 2731 ! b011 : O(1D) + HD -> H + OD 2732 !=========================================================== 2733 2734 nb_reaction_4 = nb_reaction_4 + 1 2735 2736 indice_4(nb_reaction_4) = z4spec(1.0, i_o1d, 1.0, i_hd, 1.0, i_h, 1.0, i_od) 2737 2738 2739 !=========================================================== 2740 ! b012 : O(1D) + HD -> D + OH 2741 !=========================================================== 2742 2743 nb_reaction_4 = nb_reaction_4 + 1 2744 2745 indice_4(nb_reaction_4) = z4spec(1.0, i_o1d, 1.0, i_hd, 1.0, i_d, 1.0, i_oh) 2746 2747 endif !deutchem 2748 2275 2749 !=========================================================== 2276 2750 ! c001 : O + HO2 -> OH + O2 … … 2417 2891 indice_3(nb_reaction_3) = z3spec(2.0, i_h, 1.0, i_h2, 0.0, i_dummy) 2418 2892 2893 2419 2894 !=========================================================== 2420 2895 ! d001 : NO2 + O -> NO + O2 … … 2528 3003 2529 3004 indice_4(nb_reaction_4) = z4spec(1.0, i_co, 1.0, i_o, 1.0, i_co2, 0.0, i_dummy) 3005 if(deutchem) then 3006 !=========================================================== 3007 ! f001 : OD + OH -> HDO + O 3008 !=========================================================== 3009 3010 nb_reaction_4 = nb_reaction_4 + 1 3011 3012 indice_4(nb_reaction_4) = z4spec(1.0, i_od, 1.0, i_oh, 1.0, i_hdo, 1.0, i_o) 3013 3014 !=========================================================== 3015 ! f002 : OD + H2 -> HDO + H 3016 !=========================================================== 3017 3018 nb_reaction_4 = nb_reaction_4 + 1 3019 3020 indice_4(nb_reaction_4) = z4spec(1.0, i_od, 1.0, i_h2, 1.0, i_hdo, 1.0, i_h) 3021 3022 !=========================================================== 3023 ! f003 : OD + HO2 -> HDO + O2 3024 !=========================================================== 3025 3026 nb_reaction_4 = nb_reaction_4 + 1 3027 3028 indice_4(nb_reaction_4) = z4spec(1.0, i_od, 1.0, i_ho2, 1.0, i_hdo, 1.0, i_o2) 3029 3030 !=========================================================== 3031 ! f004 : OD + H2O2 -> HDO + HO2 3032 !=========================================================== 3033 3034 nb_reaction_4 = nb_reaction_4 + 1 3035 3036 indice_4(nb_reaction_4) = z4spec(1.0, i_od, 1.0, i_h2o2, 1.0, i_hdo, 1.0, i_ho2) 3037 3038 !=========================================================== 3039 ! f005 : O + OD -> O2 + D 3040 !=========================================================== 3041 3042 nb_reaction_4 = nb_reaction_4 + 1 3043 3044 indice_4(nb_reaction_4) = z4spec(1.0, i_o, 1.0, i_od, 1.0, i_o2, 1.0, i_d) 3045 3046 !=========================================================== 3047 ! f006 : OD + H2 -> H2O + D 3048 !=========================================================== 3049 3050 nb_reaction_4 = nb_reaction_4 + 1 3051 3052 indice_4(nb_reaction_4) = z4spec(1.0, i_h2, 1.0, i_od, 1.0, i_h2o, 1.0, i_d) 3053 3054 !=========================================================== 3055 ! f007 : OD + H -> OH + D 3056 !=========================================================== 3057 3058 nb_reaction_4 = nb_reaction_4 + 1 3059 3060 indice_4(nb_reaction_4) = z4spec(1.0, i_h, 1.0, i_od, 1.0, i_oh, 1.0, i_d) 3061 3062 !=========================================================== 3063 ! f008 : CO + OD -> CO2 + D 3064 !=========================================================== 3065 3066 nb_reaction_4 = nb_reaction_4 + 1 3067 3068 indice_4(nb_reaction_4) = z4spec(1.0, i_co, 1.0, i_od, 1.0, i_co2, 1.0, i_d) 3069 3070 !=========================================================== 3071 ! f009 : O3 + D -> O2 + OD 3072 !=========================================================== 3073 3074 nb_reaction_4 = nb_reaction_4 + 1 3075 3076 indice_4(nb_reaction_4) = z4spec(1.0, i_o3, 1.0, i_d, 1.0, i_o2, 1.0, i_od) 3077 3078 !=========================================================== 3079 ! f010 : HO2 + D -> OH + OD 3080 !=========================================================== 3081 3082 nb_reaction_4 = nb_reaction_4 + 1 3083 3084 indice_4(nb_reaction_4) = z4spec(1.0, i_ho2, 1.0, i_d, 1.0, i_oh, 1.0, i_od) 3085 3086 !=========================================================== 3087 ! f011 : HO2 + D -> HDO + O 3088 !=========================================================== 3089 3090 nb_reaction_4 = nb_reaction_4 + 1 3091 3092 indice_4(nb_reaction_4) = z4spec(1.0, i_ho2, 1.0, i_d, 1.0, i_hdo, 1.0, i_o) 3093 3094 !=========================================================== 3095 ! f012 : OH + D -> H + OD 3096 !=========================================================== 3097 3098 nb_reaction_4 = nb_reaction_4 + 1 3099 3100 indice_4(nb_reaction_4) = z4spec(1.0, i_oh, 1.0, i_d, 1.0, i_h, 1.0, i_od) 3101 3102 !=========================================================== 3103 ! f013 : H + D + CO2 -> HD + CO2 3104 !=========================================================== 3105 3106 nb_reaction_4 = nb_reaction_4 + 1 3107 3108 indice_4(nb_reaction_4) = z4spec(1.0, i_h, 1.0, i_d, 1.0, i_hd, 0.0, i_dummy) 3109 3110 !=========================================================== 3111 ! f014 : D + HO2 -> HD + O2 3112 !=========================================================== 3113 3114 nb_reaction_4 = nb_reaction_4 + 1 3115 3116 indice_4(nb_reaction_4) = z4spec(1.0, i_d, 1.0, i_ho2, 1.0, i_hd, 1.0, i_o2) 3117 3118 3119 !=========================================================== 3120 ! f015 : OH + HD -> HDO + H 3121 !=========================================================== 3122 3123 nb_reaction_4 = nb_reaction_4 + 1 3124 3125 indice_4(nb_reaction_4) = z4spec(1.0, i_oh, 1.0, i_hd, 1.0, i_hdo, 1.0, i_h) 3126 3127 3128 !=========================================================== 3129 ! f016 : OH + HD -> H2O + D 3130 !=========================================================== 3131 3132 nb_reaction_4 = nb_reaction_4 + 1 3133 3134 indice_4(nb_reaction_4) = z4spec(1.0, i_oh, 1.0, i_hd, 1.0, i_h2o, 1.0, i_d) 3135 3136 3137 !=========================================================== 3138 ! f017 : D + O2 + CO2 -> DO2 + CO2 3139 !=========================================================== 3140 3141 nb_reaction_4 = nb_reaction_4 + 1 3142 3143 indice_4(nb_reaction_4) = z4spec(1.0, i_d, 1.0, i_o2, 1.0, i_do2, 0.0, i_dummy) 3144 3145 3146 !=========================================================== 3147 ! f018 : OD + O3 -> DO2 + O2 3148 !=========================================================== 3149 3150 nb_reaction_4 = nb_reaction_4 + 1 3151 3152 indice_4(nb_reaction_4) = z4spec(1.0, i_od, 1.0, i_o3, 1.0, i_do2, 1.0, i_o2) 3153 3154 3155 !=========================================================== 3156 ! f019 : D + HO2 -> DO2 + H 3157 !=========================================================== 3158 3159 nb_reaction_4 = nb_reaction_4 + 1 3160 3161 indice_4(nb_reaction_4) = z4spec(1.0, i_d, 1.0, i_ho2, 1.0, i_do2, 1.0, i_h) 3162 3163 3164 !=========================================================== 3165 ! f020 : O + DO2 -> OD + O2 3166 !=========================================================== 3167 3168 nb_reaction_4 = nb_reaction_4 + 1 3169 3170 indice_4(nb_reaction_4) = z4spec(1.0, i_o, 1.0, i_do2, 1.0, i_od, 1.0, i_o2) 3171 3172 3173 !=========================================================== 3174 ! f021 : H + DO2 -> OH + OD 3175 !=========================================================== 3176 3177 nb_reaction_4 = nb_reaction_4 + 1 3178 3179 indice_4(nb_reaction_4) = z4spec(1.0, i_h, 1.0, i_do2, 1.0, i_od, 1.0, i_oh) 3180 3181 3182 !=========================================================== 3183 ! f022 : H + DO2 -> HD + O2 3184 !=========================================================== 3185 3186 nb_reaction_4 = nb_reaction_4 + 1 3187 3188 indice_4(nb_reaction_4) = z4spec(1.0, i_h, 1.0, i_do2, 1.0, i_hd, 1.0, i_o2) 3189 3190 3191 !=========================================================== 3192 ! f023 : H + DO2 -> HDO + O 3193 !=========================================================== 3194 3195 nb_reaction_4 = nb_reaction_4 + 1 3196 3197 indice_4(nb_reaction_4) = z4spec(1.0, i_h, 1.0, i_do2, 1.0, i_hdo, 1.0, i_o) 3198 3199 3200 !=========================================================== 3201 ! f024 : H + DO2 -> HO2 + D 3202 !=========================================================== 3203 3204 nb_reaction_4 = nb_reaction_4 + 1 3205 3206 indice_4(nb_reaction_4) = z4spec(1.0, i_h, 1.0, i_do2, 1.0, i_ho2, 1.0, i_d) 3207 3208 3209 !=========================================================== 3210 ! f025 : OH + DO2 -> HDO + O2 3211 !=========================================================== 3212 3213 nb_reaction_4 = nb_reaction_4 + 1 3214 3215 indice_4(nb_reaction_4) = z4spec(1.0, i_oh, 1.0, i_do2, 1.0, i_hdo, 1.0, i_o2) 3216 3217 3218 !=========================================================== 3219 ! f026 : DO2 + O3 -> OD + O2 + O2 3220 !=========================================================== 3221 3222 nb_reaction_4 = nb_reaction_4 + 1 3223 3224 indice_4(nb_reaction_4) = z4spec(1.0, i_do2, 1.0, i_o3, 1.0, i_od, 2.0, i_o2) 3225 3226 3227 !=========================================================== 3228 ! f027 : OD + OH + CO2 -> HDO2 + CO2 3229 !=========================================================== 3230 3231 nb_reaction_4 = nb_reaction_4 + 1 3232 3233 indice_4(nb_reaction_4) = z4spec(1.0, i_od, 1.0, i_oh, 1.0, i_hdo2, 0.0, i_dummy) 3234 3235 3236 !=========================================================== 3237 ! f028 : DO2 + HO2 -> HDO2 + O2 3238 !=========================================================== 3239 3240 nb_reaction_4 = nb_reaction_4 + 1 3241 3242 indice_4(nb_reaction_4) = z4spec(1.0, i_do2, 1.0, i_ho2, 1.0, i_hdo2, 1.0, i_o2) 3243 3244 !=========================================================== 3245 ! f029 : O + HDO2 -> OD + HO2 3246 !=========================================================== 3247 3248 nb_reaction_4 = nb_reaction_4 + 1 3249 3250 indice_4(nb_reaction_4) = z4spec(1.0, i_o, 1.0, i_hdo2, 1.0, i_od, 1.0, i_ho2) 3251 3252 3253 !=========================================================== 3254 ! f030 : O + HDO2 -> OH + DO2 3255 !=========================================================== 3256 3257 nb_reaction_4 = nb_reaction_4 + 1 3258 3259 indice_4(nb_reaction_4) = z4spec(1.0, i_o, 1.0, i_hdo2, 1.0, i_oh, 1.0, i_do2) 3260 3261 3262 !=========================================================== 3263 ! f031 : OH + HDO2 -> HDO + HO2 3264 !=========================================================== 3265 3266 nb_reaction_4 = nb_reaction_4 + 1 3267 3268 indice_4(nb_reaction_4) = z4spec(1.0, i_oh, 1.0, i_hdo2, 1.0, i_hdo, 1.0, i_ho2) 3269 3270 3271 !=========================================================== 3272 ! f032 : OH + HDO2 -> H2O + DO2 3273 !=========================================================== 3274 3275 nb_reaction_4 = nb_reaction_4 + 1 3276 3277 indice_4(nb_reaction_4) = z4spec(1.0, i_oh, 1.0, i_hdo2, 1.0, i_h2o, 1.0, i_do2) 3278 3279 3280 endif !deutchem 2530 3281 2531 3282 !Only if ion chemistry … … 3091 3842 !***************************************************************** 3092 3843 3093 subroutine gcmtochim(nlayer, ionchem, nq, zycol, lswitch, nesp,& 3844 subroutine gcmtochim(nlayer, ionchem, deutchem, nq, zycol, & 3845 lswitch, nesp, & 3094 3846 i_co2, i_co, i_o, i_o1d, i_o2, i_o3, i_h, & 3095 3847 i_h2, i_oh, i_ho2, i_h2o2, i_h2o, & … … 3098 3850 i_coplus, i_cplus, i_n2plus, i_nplus, & 3099 3851 i_hplus, i_hco2plus, i_hcoplus, i_h2oplus,& 3100 i_h3oplus, i_ohplus, i_elec, dens, rm, c) 3852 i_h3oplus, i_ohplus, i_elec, i_hdo, i_od, & 3853 i_d, i_hd, i_do2, i_hdo2, dens, rm, c) 3101 3854 3102 3855 !***************************************************************** … … 3110 3863 & igcm_n2plus, igcm_nplus, igcm_hplus, & 3111 3864 & igcm_hco2plus, igcm_hcoplus, igcm_h2oplus, & 3112 & igcm_h3oplus, igcm_ohplus, igcm_elec 3865 & igcm_h3oplus, igcm_ohplus, igcm_elec, & 3866 & igcm_hdo_vap, igcm_od, igcm_d, igcm_hd, & 3867 & igcm_do2, igcm_hdo2 3113 3868 3114 3869 implicit none … … 3123 3878 integer, intent(in) :: nq ! number of tracers in the gcm 3124 3879 logical, intent(in) :: ionchem 3880 logical, intent(in) :: deutchem 3125 3881 integer :: nesp ! number of species in the chemistry 3126 3882 integer :: lswitch ! interface level between chemistries … … 3131 3887 i_co2plus, i_oplus, i_o2plus, i_noplus, i_coplus, & 3132 3888 i_cplus, i_n2plus, i_nplus, i_hplus, i_hco2plus, & 3133 i_hcoplus, i_h2oplus, i_h3oplus, i_ohplus, i_elec 3889 i_hcoplus, i_h2oplus, i_h3oplus, i_ohplus, i_elec, & 3890 i_hdo, i_od, i_d, i_hd, i_do2, i_hdo2 3134 3891 3135 3892 real :: zycol(nlayer,nq) ! volume mixing ratios in the gcm … … 3224 3981 call abort_physic("gcmtochim","missing h2o_vap tracer",1) 3225 3982 end if 3983 if(deutchem) then 3984 if (igcm_hdo_vap == 0) then 3985 write(*,*) "gcmtochim: Error; no HDO tracer !!!" 3986 call abort_physic("gcmtochim","missing hdo_vap tracer",1) 3987 end if 3988 if (igcm_od == 0) then 3989 write(*,*) "gcmtochim: Error; no OD tracer !!!" 3990 call abort_physic("gcmtochim","missing od tracer",1) 3991 end if 3992 if (igcm_d == 0) then 3993 write(*,*) "gcmtochim: Error; no D tracer !!!" 3994 call abort_physic("gcmtochim","missing d tracer",1) 3995 end if 3996 if (igcm_hd == 0) then 3997 write(*,*) "gcmtochim: Error; no HD tracer !!!" 3998 call abort_physic("gcmtochim","missing hd tracer",1) 3999 end if 4000 if (igcm_do2 == 0) then 4001 write(*,*) "gcmtochim: Error; no DO2 tracer !!!" 4002 call abort_physic("gcmtochim","missing do2 tracer",1) 4003 end if 4004 if (igcm_hdo2 == 0) then 4005 write(*,*) "gcmtochim: Error; no HDO2 tracer !!!" 4006 call abort_physic("gcmtochim","missing hdo2 tracer",1) 4007 end if 4008 endif !deutchem 3226 4009 if (ionchem) then 3227 4010 if (igcm_co2plus == 0) then … … 3311 4094 rm(l,i_no2) = zycol(l, igcm_no2) 3312 4095 rm(l,i_n2) = zycol(l, igcm_n2) 3313 end do 4096 enddo 4097 4098 if (deutchem) then 4099 do l = 1,nlayer 4100 rm(l,i_hdo) = zycol(l, igcm_hdo_vap) 4101 rm(l,i_od) = zycol(l, igcm_od) 4102 rm(l,i_d) = zycol(l, igcm_d) 4103 rm(l,i_hd) = zycol(l, igcm_hd) 4104 rm(l,i_do2) = zycol(l, igcm_do2) 4105 rm(l,i_hdo2) = zycol(l, igcm_hdo2) 4106 end do 4107 endif 3314 4108 3315 4109 if (ionchem) then … … 3351 4145 !***************************************************************** 3352 4146 3353 subroutine chimtogcm(nlayer, ionchem, nq, zycol, lswitch, nesp, & 4147 subroutine chimtogcm(nlayer, ionchem, deutchem, nq, zycol, & 4148 lswitch, nesp, & 3354 4149 i_co2, i_co, i_o, i_o1d, i_o2, i_o3, i_h, & 3355 4150 i_h2, i_oh, i_ho2, i_h2o2, i_h2o, & … … 3358 4153 i_coplus, i_cplus, i_n2plus, i_nplus, & 3359 4154 i_hplus, i_hco2plus, i_hcoplus, i_h2oplus, & 3360 i_h3oplus, i_ohplus, i_elec, dens, c) 4155 i_h3oplus, i_ohplus, i_elec, i_hdo, i_od, & 4156 i_d, i_hd, i_do2, i_hdo2, dens, c) 3361 4157 3362 4158 !***************************************************************** … … 3370 4166 igcm_n2plus, igcm_nplus, igcm_hplus, & 3371 4167 igcm_hco2plus, igcm_hcoplus, igcm_h2oplus, & 3372 igcm_h3oplus, igcm_ohplus, igcm_elec 4168 igcm_h3oplus, igcm_ohplus, igcm_elec, & 4169 igcm_hdo_vap, igcm_od, igcm_d, igcm_hd, & 4170 igcm_do2, igcm_hdo2 3373 4171 3374 4172 implicit none … … 3383 4181 integer, intent(in) :: nq ! number of tracers in the gcm 3384 4182 logical, intent(in) :: ionchem 4183 logical, intent(in) :: deutchem 3385 4184 integer :: nesp ! number of species in the chemistry 3386 4185 integer :: lswitch ! interface level between chemistries … … 3391 4190 i_coplus, i_cplus, i_n2plus, i_nplus, & 3392 4191 i_hplus, i_hco2plus, i_hcoplus, i_h2oplus, & 3393 i_h3oplus, i_ohplus, i_elec 4192 i_h3oplus, i_ohplus, i_elec, i_hdo, i_od, i_d, & 4193 i_hd, i_do2, i_hdo2 3394 4194 3395 4195 real :: dens(nlayer) ! total number density (molecule.cm-3) … … 3430 4230 zycol(l, igcm_no2) = c(l,i_no2)/dens(l) 3431 4231 zycol(l, igcm_n2) = c(l,i_n2)/dens(l) 3432 end do 4232 enddo 4233 4234 if(deutchem) then 4235 do l=1,lswitch-1 4236 zycol(l, igcm_hdo_vap) = c(l,i_hdo)/dens(l) 4237 zycol(l, igcm_od) = c(l,i_od)/dens(l) 4238 zycol(l, igcm_d) = c(l,i_d)/dens(l) 4239 zycol(l, igcm_hd) = c(l,i_hd)/dens(l) 4240 zycol(l, igcm_do2) = c(l,i_do2)/dens(l) 4241 zycol(l, igcm_hdo2) = c(l,i_hdo2)/dens(l) 4242 end do 4243 endif !deutchem 3433 4244 3434 4245 if (ionchem) then -
trunk/LMDZ.MARS/libf/aeronomars/photolysis_mod.F90
r2433 r2461 5 5 ! photolysis 6 6 7 integer, parameter:: nphot = 13 ! number of photolysis8 ! integer, parameter :: nphot = 14! number of photolysis (with hdo)7 integer, save :: nphot = 13 ! number of photolysis 8 ! integer, parameter :: nphot = 15 ! number of photolysis (with hdo) 9 9 integer, parameter :: nabs = 10 ! number of absorbing gases 10 10 -
trunk/LMDZ.MARS/libf/aeronomars/photolysis_online.F
r2433 r2461 1 1 !============================================================================== 2 2 3 subroutine photolysis_online(nlayer, nb_phot_max, alt, press, 4 $ temp, zmmean, i_co2, i_co, i_o, i_o1d, i_o2, i_o3, i_h, 3 subroutine photolysis_online(nlayer, deutchem, nb_phot_max, 4 $ alt, press, temp, zmmean, 5 $ i_co2, i_co, i_o, i_o1d, i_o2, i_o3, i_h, 5 6 $ i_h2, i_oh, i_ho2, i_h2o2, i_h2o, 6 7 $ i_n, i_n2d, i_no, i_no2, i_n2, nesp, rm, … … 13 14 ! input 14 15 16 logical, intent(in) :: deutchem 15 17 integer, intent(in) :: nesp ! total number of chemical species 16 18 integer, intent(in) :: nlayer … … 60 62 61 63 integer :: j_o2_o, j_o2_o1d, j_co2_o, j_co2_o1d, j_o3_o1d, j_o3_o, 62 $ j_h2o, j_h2o2, j_ho2, j_h2, j_no, j_no2, j_n2, j_hdo 64 $ j_h2o, j_h2o2, j_ho2, j_h2, j_no, j_no2, j_n2, 65 $ j_hdo_od, j_hdo_d 63 66 64 67 integer :: a_o2, a_co2, a_o3, a_h2o, a_h2o2, a_ho2, a_h2, a_no, … … 96 99 j_no2 = 12 ! no2 + hv -> no + o 97 100 j_n2 = 13 ! n2 + hv -> n + n 98 j_hdo = 14 ! hdo + hv -> products 101 j_hdo_od = 14 ! hdo + hv -> od + h 102 j_hdo_d = 15 ! hdo + hv -> d + oh 99 103 100 104 !==== air column increments and rayleigh optical depth … … 168 172 sj(ilay,iw,j_no) = xsno(iw)*yieldno(iw) ! no 169 173 sj(ilay,iw,j_n2) = xsn2(iw)*yieldn2(iw) ! n2 170 ! sj(ilay,iw,j_hdo) = xshdo(iw) ! hdo171 174 end do 172 175 end do 176 177 !HDO cross section 178 if (deutchem) then 179 do ilay=1,nlayer 180 do iw=1,nw-1 181 !Two chanels for HDO, with same cross section 182 sj(ilay,iw,j_hdo_od) = 0.5*xshdo(iw) ! hdo -> od + h 183 sj(ilay,iw,j_hdo_d) = 0.5*xshdo(iw) ! hdo -> d + oh 184 enddo 185 enddo 186 endif 173 187 174 188 !==== set aerosol properties and optical depth
Note: See TracChangeset
for help on using the changeset viewer.
