Changeset 2170 for trunk/LMDZ.MARS/libf/aeronomars/calchim_mod.F90
- Timestamp:
- Oct 14, 2019, 11:00:51 PM (5 years ago)
- File:
-
- 1 edited
-
trunk/LMDZ.MARS/libf/aeronomars/calchim_mod.F90 (modified) (12 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/LMDZ.MARS/libf/aeronomars/calchim_mod.F90
r2164 r2170 26 26 use conc_mod, only: mmean ! mean molecular mass of the atmosphere 27 27 use comcstfi_h, only: pi 28 use photolysis_mod, only: jonline, init_photolysis28 use photolysis_mod, only: init_photolysis, nphot 29 29 use iono_h, only: temp_elect 30 30 … … 38 38 ! Prepare the call for the photochemical module, and send back the 39 39 ! tendencies from photochemistry in the chemical species mass mixing ratios 40 !41 ! Author: Sebastien Lebonnois (08/11/2002)42 ! -------43 ! update 12/06/2003 for water ice clouds and compatibility with dust44 ! update 07/2003 for coupling with thermosphere (Monica Angelats-i-Coll)45 ! update 03/05/2005 cosmetic changes (Franck Lefevre)46 ! update sept. 2008 identify tracers by their names (Ehouarn Millour)47 ! update 05/12/2011 synchronize with latest version of chemistry (Franck Lefevre)48 ! update 16/03/2012 optimization (Franck Lefevre)49 ! update 11/12/2013 optimization (Franck Lefevre)50 40 ! 51 41 ! Arguments: … … 155 145 integer :: ig_vl1 156 146 147 integer :: nb_reaction_3_max ! number of quadratic reactions 148 integer :: nb_reaction_4_max ! number of bimolecular reactions 149 integer :: nquench ! number of quenching + heterogeneous reactions 150 integer :: nphotion ! number of photoionizations 151 integer :: nb_phot_max ! total number of photolysis+photoionizations+quenching reactions 152 153 157 154 real :: latvl1, lonvl1 158 155 real :: zq(ngrid,nlayer,nq) ! pq+pdq*ptimestep before chemistry … … 164 161 real :: kb ! boltzmann constant 165 162 166 logical,save :: firstcall = .true. 167 logical,save :: depos = .false. ! switch for dry deposition 163 logical, save :: firstcall = .true. 164 logical, save :: depos ! switch for dry deposition 165 logical, save :: ionchem ! switch for ionospheric chemistry 166 logical, save :: jonline ! switch for online photodissociation rates or lookup table 167 logical, save :: unichim ! only one unified chemical scheme at all 168 ! layers (default), or upper atmospheric 169 ! scheme in the thermosphere 168 170 169 171 ! for each column of atmosphere: … … 181 183 real :: jo3(nlayer) ! Photodissociation rate O3->O1D (s-1) 182 184 real :: jh2o(nlayer) ! Photodissociation rate H2O->H+OH (s-1) 183 real :: em_no(nlayer) ! 184 real :: em_o2(nlayer) ! 185 186 integer :: iter(nlayer) ! 187 ! 185 real :: em_no(nlayer) ! NO nightglow emission rate 186 real :: em_o2(nlayer) ! O2 nightglow emission rate 187 188 integer :: iter(nlayer) ! Number of chemical iterations 189 ! within one physical timestep 188 190 189 191 ! for output: 190 192 191 193 logical :: output ! to issue calls to writediagfi and stats 192 parameter (output = .true.)193 194 real :: jo3_3d(ngrid,nlayer) ! Photodissociation rate O3->O1D (s-1) 194 195 real :: jh2o_3d(ngrid,nlayer) ! Photodissociation rate H2O->H+OH (s-1) … … 198 199 ! within one physical timestep 199 200 200 logical :: unichim ! only one, unified chemical scheme at all 201 ! layers (default), or upper atmospheric 202 ! scheme in the thermosphere? 201 !======================================================================= 202 ! main dashboard for the chemistry 203 !======================================================================= 204 205 unichim = .true. ! true : unified chemistry ! false : separate models in lower and upper atmosphere 206 jonline = .true. ! true : on-line calculation of photodissociation rates ! false : lookup table 207 ionchem = .false. ! switch for ionospheric chemistry 208 depos = .false. ! switch for dry deposition 209 output = .false. ! true : triggers writing of specific outputs (photolysis rates, emission rates, etc) 203 210 204 211 !======================================================================= 205 212 ! initialization of the chemistry (first call only) 206 213 !======================================================================= 207 208 !call only unified chemistry (default), or also upper atmospheric model209 !unichim = .true. unified chemistry210 !unichim = .false. 2 different models211 unichim = .true.212 214 213 215 if (firstcall) then … … 388 390 i_o2plus = igcm_o2plus 389 391 if (i_o2plus == 0) then 390 write(*,*) "calchim: Error, no O2+ tracer !!!"391 stop392 write(*,*) "calchim: no O2+ tracer " 393 write(*,*) "Only neutral chemistry" 392 394 else 393 395 nbq = nbq + 1 394 396 niq(nbq) = i_o2plus 397 ionchem = .true. 398 write(*,*) "calchim: O2+ tracer found in traceur.def" 399 write(*,*) "Ion chemistry included" 395 400 end if 396 401 i_co2plus = igcm_co2plus 397 if (i_co2plus == 0) then 398 write(*,*) "calchim: Error, no CO2+ tracer !!!" 399 stop 400 else 401 nbq = nbq + 1 402 niq(nbq) = i_co2plus 403 end if 402 if(ionchem) then 403 if (i_co2plus == 0) then 404 write(*,*) "calchim: Error, no CO2+ tracer !!!" 405 write(*,*) "CO2+ is needed if O2+ is in traceur.def" 406 stop 407 else 408 nbq = nbq + 1 409 niq(nbq) = i_co2plus 410 end if 411 else 412 if (i_co2plus /= 0) then 413 write(*,*) "calchim: Error: CO2+ is present, but O2+ is not!!!" 414 write(*,*) "Both must be in traceur.def if ion chemistry wanted" 415 stop 416 endif 417 endif 404 418 i_oplus=igcm_oplus 405 if (i_oplus == 0) then 406 write(*,*) "calchim: Error, no O+ tracer !!!" 407 stop 408 else 409 nbq = nbq + 1 410 niq(nbq) = i_oplus 411 end if 419 if(ionchem) then 420 if (i_oplus == 0) then 421 write(*,*) "calchim: Error, no O+ tracer !!!" 422 write(*,*) "O+ is needed if O2+ is in traceur.def" 423 stop 424 else 425 nbq = nbq + 1 426 niq(nbq) = i_oplus 427 end if 428 else 429 if (i_oplus /= 0) then 430 write(*,*) "calchim: Error: O+ is present, but O2+ is not!!!" 431 write(*,*) "Both must be in traceur.def if ion chemistry wanted" 432 stop 433 endif 434 endif 412 435 i_noplus=igcm_noplus 413 if (i_noplus == 0) then 414 write(*,*) "calchim: Error, no NO+ tracer !!!" 415 stop 416 else 417 nbq = nbq + 1 418 niq(nbq) = i_noplus 419 end if 436 if(ionchem) then 437 if (i_noplus == 0) then 438 write(*,*) "calchim: Error, no NO+ tracer !!!" 439 write(*,*) "NO+ is needed if O2+ is in traceur.def" 440 stop 441 else 442 nbq = nbq + 1 443 niq(nbq) = i_noplus 444 end if 445 else 446 if (i_noplus /= 0) then 447 write(*,*) "calchim: Error: NO+ is present, but O2+ is not!!!" 448 write(*,*) "Both must be in traceur.def if ion chemistry wanted" 449 endif 450 endif 420 451 i_coplus=igcm_coplus 421 if (i_coplus == 0) then 422 write(*,*) "calchim: Error, no CO+ tracer !!!" 423 stop 424 else 425 nbq = nbq + 1 426 niq(nbq) = i_coplus 427 end if 452 if(ionchem) then 453 if (i_coplus == 0) then 454 write(*,*) "calchim: Error, no CO+ tracer !!!" 455 write(*,*) "CO+ is needed if O2+ is in traceur.def" 456 stop 457 else 458 nbq = nbq + 1 459 niq(nbq) = i_coplus 460 end if 461 else 462 if (i_coplus /= 0) then 463 write(*,*) "calchim: Error: CO+ is present, but O2+ is not!!!" 464 write(*,*) "Both must be in traceur.def if ion chemistry wanted" 465 endif 466 endif 428 467 i_cplus=igcm_cplus 429 if (i_cplus == 0) then 430 write(*,*) "calchim: Error, no C+ tracer !!!" 431 stop 432 else 433 nbq = nbq + 1 434 niq(nbq) = i_cplus 435 end if 468 if(ionchem) then 469 if (i_cplus == 0) then 470 write(*,*) "calchim: Error, no C+ tracer !!!" 471 write(*,*) "C+ is needed if O2+ is in traceur.def" 472 stop 473 else 474 nbq = nbq + 1 475 niq(nbq) = i_cplus 476 end if 477 else 478 if (i_cplus /= 0) then 479 write(*,*) "calchim: Error: C+ is present, but O2+ is not!!!" 480 write(*,*) "Both must be in traceur.def if ion chemistry wanted" 481 endif 482 endif 436 483 i_n2plus=igcm_n2plus 437 if (i_n2plus == 0) then 438 write(*,*) "calchim: Error, no N2+ tracer !!!" 439 stop 440 else 441 nbq = nbq + 1 442 niq(nbq) = i_n2plus 443 end if 484 if(ionchem) then 485 if (i_n2plus == 0) then 486 write(*,*) "calchim: Error, no N2+ tracer !!!" 487 write(*,*) "N2+ is needed if O2+ is in traceur.def" 488 stop 489 else 490 nbq = nbq + 1 491 niq(nbq) = i_n2plus 492 end if 493 else 494 if (i_n2plus /= 0) then 495 write(*,*) "calchim: Error: N2+ is present, but O2+ is not!!!" 496 write(*,*) "Both must be in traceur.def if ion chemistry wanted" 497 endif 498 endif 444 499 i_nplus=igcm_nplus 445 if (i_nplus == 0) then 446 write(*,*) "calchim: Error, no N+ tracer !!!" 447 stop 448 else 449 nbq = nbq + 1 450 niq(nbq) = i_nplus 451 end if 500 if(ionchem) then 501 if (i_nplus == 0) then 502 write(*,*) "calchim: Error, no N+ tracer !!!" 503 write(*,*) "N+ is needed if O2+ is in traceur.def" 504 stop 505 else 506 nbq = nbq + 1 507 niq(nbq) = i_nplus 508 end if 509 else 510 if (i_nplus /= 0) then 511 write(*,*) "calchim: Error: N+ is present, but O2+ is not!!!" 512 write(*,*) "Both must be in traceur.def if ion chemistry wanted" 513 endif 514 endif 452 515 i_hplus=igcm_hplus 453 if (i_hplus == 0) then 454 write(*,*) "calchim: Error, no H+ tracer !!!" 455 stop 456 else 457 nbq = nbq + 1 458 niq(nbq) = i_hplus 459 end if 516 if(ionchem) then 517 if (i_hplus == 0) then 518 write(*,*) "calchim: Error, no H+ tracer !!!" 519 write(*,*) "H+ is needed if O2+ is in traceur.def" 520 stop 521 else 522 nbq = nbq + 1 523 niq(nbq) = i_hplus 524 end if 525 else 526 if (i_hplus /= 0) then 527 write(*,*) "calchim: Error: H+ is present, but O2+ is not!!!" 528 write(*,*) "Both must be in traceur.def if ion chemistry wanted" 529 endif 530 endif 460 531 i_hco2plus=igcm_hco2plus 461 if (i_hco2plus == 0) then 462 write(*,*) "calchim: Error, no HCO2+ tracer !!!" 463 stop 464 else 465 nbq = nbq + 1 466 niq(nbq) = i_hco2plus 467 end if 532 if(ionchem) then 533 if (i_hco2plus == 0) then 534 write(*,*) "calchim: Error, no HCO2+ tracer !!!" 535 write(*,*) "HCO2+ is needed if O2+ is in traceur.def" 536 stop 537 else 538 nbq = nbq + 1 539 niq(nbq) = i_hco2plus 540 end if 541 else 542 if (i_hco2plus /= 0) then 543 write(*,*) "calchim: Error: HCO2+ is present, but O2+ is not!!!" 544 write(*,*) "Both must be in traceur.def if ion chemistry wanted" 545 endif 546 endif 468 547 i_elec = igcm_elec 469 if (i_elec == 0) then 470 write(*,*) "calchim: Error, no e- tracer !!!" 471 stop 472 else 473 nbq = nbq + 1 474 niq(nbq) = i_elec 475 end if 476 548 if(ionchem) then 549 if (i_elec == 0) then 550 write(*,*) "calchim: Error, no e- tracer !!!" 551 write(*,*) "e- is needed if O2+ is in traceur.def" 552 stop 553 else 554 nbq = nbq + 1 555 niq(nbq) = i_elec 556 end if 557 else 558 if (i_elec /= 0) then 559 write(*,*) "calchim: Error: e- is present, but O2+ is not!!!" 560 write(*,*) "Both must be in traceur.def if ion chemistry wanted" 561 endif 562 endif 477 563 write(*,*) 'calchim: found nbq = ',nbq,' tracers' 478 564 … … 523 609 524 610 ! search for switch index between regions 525 if(unichim) then 526 lswitch=nlayer+1 611 612 if (unichim) then 613 lswitch = nlayer + 1 527 614 else 528 615 if (foundswitch == 0 .and. pplay(ig,l) < 1.e-1) then … … 540 627 !======================================================================= 541 628 542 ! chemistry: only one scheme at all layers543 544 629 if (photochem) then 545 call photochemistry(nlayer,nq, & 546 ig,lswitch,zycol,szacol,ptimestep, & 547 zpress,zlocal,ztemp,ztelec,zdens,zmmean, & 548 dist_sol,zday,surfdust1d,surfice1d, & 630 ! set number of reactions, depending on ion chemistry or not 631 if (ionchem) then 632 nb_reaction_4_max = 60 ! set number of bimolecular reactions 633 nb_reaction_3_max = 6 ! set number of quadratic reactions 634 nquench = 9 ! set number of quenching + heterogeneous reactions 635 nphotion = 18 ! set number of photoionizations 636 else 637 nb_reaction_4_max = 31 ! set number of bimolecular reactions 638 nb_reaction_3_max = 6 ! set number of quadratic reactions 639 nquench = 9 ! set number of quenching + heterogeneous reactions 640 nphotion = 0 ! set number of photoionizations 641 end if 642 643 ! nb_phot_max is the total number of processes that are treated 644 ! numerically as a photolysis: 645 646 nb_phot_max = nphot + nphotion + nquench 647 648 ! call main photochemistry routine 649 650 call photochemistry(nlayer,nq,nbq,ionchem,nb_reaction_3_max, & 651 nb_reaction_4_max, nb_phot_max, nphotion, & 652 jonline, ig,lswitch,zycol,szacol,ptimestep, & 653 zpress,zlocal,ztemp,ztelec,zdens,zmmean, & 654 dist_sol,zday,surfdust1d,surfice1d, & 549 655 jo3,jh2o,taucol,iter) 550 656 … … 556 662 iter_3d(ig,l) = iter(l) 557 663 end do 664 558 665 ! condensation of h2o2 559 666 … … 562 669 ztemp,zycol,dqcloud,dqscloud) 563 670 564 if(.not.unichim) then 565 chemthermod=3 !C/O/H/N/ions 671 ! case of separate photochemical model in the thermosphere 672 673 if (.not.unichim) then 674 chemthermod = 3 !C/O/H/N/ions 566 675 call chemthermos(ig,nlayer,lswitch,chemthermod,zycol,ztemp,& 567 676 zdens,zpress,zlocal,szacol,ptimestep,zday,& 568 677 em_no,em_o2) 569 do l=1,nlayer 570 emission_no(ig,l)=em_no(l) 571 emission_o2(ig,l)=em_o2(l) 572 enddo 573 end if 574 end if 678 do l = 1,nlayer 679 emission_no(ig,l) = em_no(l) 680 emission_o2(ig,l) = em_o2(l) 681 end do 682 end if 683 684 end if ! photochem 575 685 576 686 ! dry deposition … … 627 737 call writediagfi(ngrid,'iter','iterations', & 628 738 ' ',3,iter_3d(1,1)) 629 if(.not.unichim) then 739 740 if (.not. unichim) then 630 741 call writediagfi(ngrid,'emission_no', & 631 742 'NO nightglow emission rate','cm-3 s-1',3,emission_no)
Note: See TracChangeset
for help on using the changeset viewer.
