source: trunk/LMDZ.MARS/libf/aeronomars/calchim_mod.F90 @ 3165

Last change on this file since 3165 was 3012, checked in by emillour, 2 years ago

Mars PCM:
Code cleanup concerning chemistry. Turn chimiedata.h into module
chemistrydata.F90 and integrate read_phototable.F in it.
Also fix an OpenMP issue with read_phototable; different OpenMP threads
should not simultaneously open/read a file. Let the master do the job
and then broadcast the result to all cores.
While at it, also turned nltecool.F and nlte_tcool.F into modules.
EM

File size: 38.5 KB
RevLine 
[2162]1MODULE calchim_mod
2
3  IMPLICIT NONE
4
[2164]5  INTEGER, SAVE :: ichemistry ! compute chemistry every ichemistry physics step
[2162]6  REAL,SAVE,ALLOCATABLE :: zdqchim(:,:,:) ! Tendancy on pq due to photochemistry
7  REAL,SAVE,ALLOCATABLE :: zdqschim(:,:) ! Tendancy on qsurf due to photochemistry
8
[2578]9!$OMP THREADPRIVATE(ichemistry,zdqchim,zdqschim)
10
[2162]11  CONTAINS
12
[2158]13      subroutine calchim(ngrid,nlayer,nq,                           &
[1495]14                         ptimestep,pplay,pplev,pt,pdt,dist_sol,mu0,         &
15                         zzlev,zzlay,zday,pq,pdq,dqchim,dqschim,dqcloud,    &
[2031]16                         dqscloud,tau,co2ice,                               &
[1495]17                         pu,pdu,pv,pdv,surfdust,surfice)
18
19      use tracer_mod, only: igcm_co2, igcm_co, igcm_o, igcm_o1d, igcm_o2, &
20                            igcm_o3, igcm_h, igcm_h2, igcm_oh, igcm_ho2,  &
21                            igcm_h2o2, igcm_ch4, igcm_n2, igcm_h2o_vap,   &
22                            igcm_no, igcm_n, igcm_no2, igcm_n2d,          &
23                            igcm_o2plus, igcm_co2plus, igcm_oplus,        &
24                            igcm_coplus, igcm_cplus, igcm_nplus,          &
25                            igcm_noplus, igcm_n2plus, igcm_hplus,         &
[2302]26                            igcm_hco2plus, igcm_hcoplus, igcm_h2oplus,    &
[2461]27                            igcm_h3oplus, igcm_ohplus, igcm_elec,         &
28                            igcm_hdo_vap, igcm_od, igcm_d, igcm_hd,       &
29                            igcm_do2, igcm_hdo2, mmol
[1495]30
31      use conc_mod, only: mmean ! mean molecular mass of the atmosphere
[2162]32      use comcstfi_h, only: pi
[3012]33      use chemistrydata_mod, only: read_phototable
[2170]34      use photolysis_mod, only: init_photolysis, nphot
[2158]35      use iono_h, only: temp_elect
[2563]36      use wstats_mod, only: wstats
[1495]37
38      implicit none
39
40!=======================================================================
41!
42!   subject:
43!   --------
44!
45!  Prepare the call for the photochemical module, and send back the
46!  tendencies from photochemistry in the chemical species mass mixing ratios
47!
48!   Arguments:
49!   ----------
50!
51!  Input:
52!
[2031]53!    ptimestep              timestep (s)
[1495]54!    pplay(ngrid,nlayer)    Pressure at the middle of the layers (Pa)
55!    pplev(ngrid,nlayer+1)  Intermediate pressure levels (Pa)
56!    pt(ngrid,nlayer)       Temperature (K)
57!    pdt(ngrid,nlayer)      Temperature tendency (K)
58!    pu(ngrid,nlayer)       u component of the wind (ms-1)
59!    pdu(ngrid,nlayer)      u component tendency (K)
60!    pv(ngrid,nlayer)       v component of the wind (ms-1)
61!    pdv(ngrid,nlayer)      v component tendency (K)
[2031]62!    dist_sol               distance of the sun (AU)
63!    mu0(ngrid)             cos of solar zenith angle (=1 when sun at zenith)
64!    pq(ngrid,nlayer,nq)    advected fields, ie chemical species here
65!    pdq(ngrid,nlayer,nq)   previous tendencies on pq
66!    tau(ngrid)             dust optical depth
67!    co2ice(ngrid)          co2 ice surface layer (kg.m-2)
[1495]68!    surfdust(ngrid,nlayer) dust surface area (m2/m3)
69!    surfice(ngrid,nlayer)  ice surface area (m2/m3)
70!
71!  Output:
72!
[2031]73!    dqchim(ngrid,nlayer,nq) tendencies on pq due to chemistry
74!    dqschim(ngrid,nq)       tendencies on qsurf
[1495]75!
76!=======================================================================
77
[2162]78include "callkeys.h"
[1495]79
80!     input:
81
82      integer,intent(in) :: ngrid    ! number of atmospheric columns
83      integer,intent(in) :: nlayer   ! number of atmospheric layers
84      integer,intent(in) :: nq       ! number of tracers
85      real :: ptimestep
86      real :: pplay(ngrid,nlayer)    ! pressure at the middle of the layers
87      real :: zzlay(ngrid,nlayer)    ! pressure at the middle of the layers
88      real :: pplev(ngrid,nlayer+1)  ! intermediate pressure levels
89      real :: zzlev(ngrid,nlayer+1)  ! altitude at layer boundaries
90      real :: pt(ngrid,nlayer)       ! temperature
91      real :: pdt(ngrid,nlayer)      ! temperature tendency
92      real :: pu(ngrid,nlayer)       ! u component of the wind (m.s-1)
93      real :: pdu(ngrid,nlayer)      ! u component tendency
94      real :: pv(ngrid,nlayer)       ! v component of the wind (m.s-1)
95      real :: pdv(ngrid,nlayer)      ! v component tendency
96      real :: dist_sol               ! distance of the sun (AU)
97      real :: mu0(ngrid)             ! cos of solar zenith angle (=1 when sun at zenith)
98      real :: pq(ngrid,nlayer,nq)    ! tracers mass mixing ratio
99      real :: pdq(ngrid,nlayer,nq)   ! previous tendencies
100      real :: zday                   ! date (time since Ls=0, in martian days)
[2031]101      real :: tau(ngrid)             ! dust optical depth
[1495]102      real :: co2ice(ngrid)          ! co2 ice surface layer (kg.m-2)
103      real :: surfdust(ngrid,nlayer) ! dust surface area (m2/m3)
[2031]104      real :: surfice(ngrid,nlayer)  ! ice surface area (m2/m3)
[1495]105
106!     output:
107
108      real :: dqchim(ngrid,nlayer,nq)   ! tendencies on pq due to chemistry
109      real :: dqschim(ngrid,nq)         ! tendencies on qsurf
110      real :: dqcloud(ngrid,nlayer,nq)  ! tendencies on pq due to condensation
111      real :: dqscloud(ngrid,nq)        ! tendencies on qsurf
112
113!     local variables:
114
115      integer,save :: nbq                   ! number of tracers used in the chemistry
116      integer,allocatable,save :: niq(:)    ! array storing the indexes of the tracers
117      integer :: iloc(1)                    ! index of major species
118      integer :: ig,l,i,iq,iqmax
119      integer :: foundswitch, lswitch
120      integer,save :: chemthermod
121
[2613]122!$OMP THREADPRIVATE(nbq,niq,chemthermod)
123
[1495]124      integer,save :: i_co2  = 0
125      integer,save :: i_co   = 0
126      integer,save :: i_o    = 0
127      integer,save :: i_o1d  = 0
128      integer,save :: i_o2   = 0
129      integer,save :: i_o3   = 0
130      integer,save :: i_h    = 0
131      integer,save :: i_h2   = 0
132      integer,save :: i_oh   = 0
133      integer,save :: i_ho2  = 0
134      integer,save :: i_h2o2 = 0
135      integer,save :: i_ch4  = 0
136      integer,save :: i_n2   = 0
137      integer,save :: i_h2o  = 0
138      integer,save :: i_n    = 0
139      integer,save :: i_no   = 0
140      integer,save :: i_no2  = 0
141      integer,save :: i_n2d  = 0
142      integer,save :: i_co2plus=0
143      integer,save :: i_oplus=0
144      integer,save :: i_o2plus=0
145      integer,save :: i_coplus=0
146      integer,save :: i_cplus=0
147      integer,save :: i_nplus=0
148      integer,save :: i_noplus=0
149      integer,save :: i_n2plus=0
150      integer,save :: i_hplus=0
151      integer,save :: i_hco2plus=0
[2284]152      integer,save :: i_hcoplus=0
[2302]153      integer,save :: i_h2oplus=0
[2321]154      integer,save :: i_h3oplus=0
155      integer,save :: i_ohplus=0
[1495]156      integer,save :: i_elec=0
[2461]157      integer,save :: i_hdo=0
158      integer,save :: i_od=0
159      integer,save :: i_d=0
160      integer,save :: i_hd=0
161      integer,save :: i_do2=0
162      integer,save :: i_hdo2=0
[1495]163
[2613]164!$OMP THREADPRIVATE(i_co2,i_co,i_o,i_o1d,i_o2,i_o3,i_h,i_h2,i_oh,i_ho2,i_h2o2,i_ch4)
165!$OMP THREADPRIVATE(i_n2,i_h2o,i_n,i_no,i_no2,i_n2d,i_co2plus,i_oplus,i_o2plus,i_coplus,i_cplus,i_nplus)
166!$OMP THREADPRIVATE(i_noplus,i_n2plus,i_hplus,i_hco2plus,i_hcoplus,i_h2oplus,i_h3oplus,i_ohplus,i_elec,i_hdo,i_od,i_d,i_hd,i_do2,i_hdo2)
167
[1495]168      integer :: ig_vl1
169
[2170]170      integer :: nb_reaction_3_max     ! number of quadratic reactions
171      integer :: nb_reaction_4_max     ! number of bimolecular reactions
172      integer :: nquench               ! number of quenching + heterogeneous reactions
173      integer :: nphotion              ! number of photoionizations
[2461]174      integer :: nb_reaction_4_ion     ! quadratic reactions for ionosphere
175      integer :: nb_reaction_4_deut    ! quadratic reactions for deuterium chem
[2170]176      integer :: nb_phot_max           ! total number of photolysis+photoionizations+quenching reactions
177
178
[1495]179      real    :: latvl1, lonvl1
180      real    :: zq(ngrid,nlayer,nq)   ! pq+pdq*ptimestep before chemistry
181                                       ! new mole fraction after
182      real    :: zt(ngrid,nlayer)      ! temperature
183      real    :: zu(ngrid,nlayer)      ! u component of the wind
184      real    :: zv(ngrid,nlayer)      ! v component of the wind
[2027]185      real    :: taucol                ! dust optical depth at the surface
[2030]186      real    :: kb                    ! boltzmann constant
[1495]187
[2170]188      logical, save :: firstcall = .true.
189      logical, save :: depos       ! switch for dry deposition
190      logical, save :: ionchem     ! switch for ionospheric chemistry
[2461]191      logical, save :: deutchem    ! switch for deuterium chemistry
[2170]192      logical, save :: jonline     ! switch for online photodissociation rates or lookup table
193      logical, save :: unichim     ! only one unified chemical scheme at all
194                                   ! layers (default), or upper atmospheric
195                                   ! scheme in the thermosphere
[1495]196
[2613]197!$OMP THREADPRIVATE(firstcall,depos,ionchem,deutchem,jonline,unichim)
198
[1495]199!     for each column of atmosphere:
200
201      real :: zpress(nlayer)       ! Pressure (mbar)
202      real :: zdens(nlayer)        ! Density  (cm-3)
203      real :: ztemp(nlayer)        ! Temperature (K)
[2158]204      real :: ztelec(nlayer)       ! Electronic temperature (K)
[1495]205      real :: zlocal(nlayer)       ! Altitude (km)
206      real :: zycol(nlayer,nq)     ! Composition (mole fractions)
207      real :: zmmean(nlayer)       ! Mean molecular mass (g.mole-1)
208      real :: szacol               ! Solar zenith angle
209      real :: surfice1d(nlayer)    ! Ice surface area (cm2/cm3)
210      real :: surfdust1d(nlayer)   ! Dust surface area (cm2/cm3)
211      real :: jo3(nlayer)          ! Photodissociation rate O3->O1D (s-1)
[2030]212      real :: jh2o(nlayer)         ! Photodissociation rate H2O->H+OH (s-1)
[2170]213      real :: em_no(nlayer)        ! NO nightglow emission rate
214      real :: em_o2(nlayer)        ! O2 nightglow emission rate     
[1495]215
[2170]216      integer :: iter(nlayer)      ! Number of chemical iterations
217                                   ! within one physical timestep
[1495]218
219!     for output:
220
[2213]221      logical,save :: output        ! to issue calls to writediagfi and stats
[1495]222      real :: jo3_3d(ngrid,nlayer)  ! Photodissociation rate O3->O1D (s-1)
[2030]223      real :: jh2o_3d(ngrid,nlayer)  ! Photodissociation rate H2O->H+OH (s-1)
[1888]224      real :: emission_no(ngrid,nlayer) !NO emission rate
225      real :: emission_o2(ngrid,nlayer) !O2 emission rate
[1495]226      real :: iter_3d(ngrid,nlayer) ! Number of chemical iterations
227                                    !  within one physical timestep
228
[2613]229!$OMP THREADPRIVATE(output)
230
[2170]231!=======================================================================
[2213]232!     initialization of the chemistry (first call only)
[2170]233!=======================================================================
[2158]234
[2213]235      if (firstcall) then
[2170]236
[1495]237!=======================================================================
[2213]238!     main dashboard for the chemistry
[1495]239!=======================================================================
240
[2213]241        unichim = .true.     ! true : unified chemistry ! false : separate models in lower and upper atmosphere
242        jonline = .true.     ! true : on-line calculation of photodissociation rates ! false : lookup table
243        ionchem = .false.    ! switch for ionospheric chemistry
[2461]244        deutchem= .false.    ! switch for deuterium chemistry
[2213]245        depos   = .false.    ! switch for dry deposition
[2461]246        output  = .true.    ! true : triggers writing of specific outputs (photolysis rates, emission rates, etc)
[1495]247
[2461]248!         if (photochem) then
249!            if (jonline) then
250!               print*,'calchim: Read UV absorption cross-sections'
251!               call init_photolysis     ! on-line photolysis
252!            else
253!               print*,'calchim: Read photolysis lookup table'
254!               call read_phototable     ! off-line photolysis
255!            end if
256!         end if
[2158]257
[2461]258!         if(.not.unichim) then
[2158]259            !Read reaction rates from external file if the upper atmospheric
260            !chemistry is called
[2461]261!            call chemthermos_readini
262!         endif
[2158]263
[1495]264         ! find index of chemical tracers to use
265         allocate(niq(nq))
266         ! Listed here are all tracers that can go into photochemistry
267         nbq = 0 ! to count number of tracers
[2158]268         ! Species ALWAYS present if photochem=.T.
[1495]269         i_co2 = igcm_co2
270         if (i_co2 == 0) then
271            write(*,*) "calchim: Error; no CO2 tracer !!!"
[2302]272            call abort_physic("calchim","missing co2 tracer",1)
[1495]273         else
274            nbq = nbq + 1
275            niq(nbq) = i_co2
276         end if
277         i_co = igcm_co
278         if (i_co == 0) then
279            write(*,*) "calchim: Error; no CO tracer !!!"
[2302]280            call abort_physic("calchim","missing co tracer",1)
[1495]281         else
282            nbq = nbq + 1
283            niq(nbq) = i_co
284         end if
285         i_o = igcm_o
286         if (i_o == 0) then
287            write(*,*) "calchim: Error; no O tracer !!!"
[2302]288            call abort_physic("calchim","missing o tracer",1)
[1495]289         else
290            nbq = nbq + 1
291            niq(nbq) = i_o
292         end if
293         i_o1d = igcm_o1d
294         if (i_o1d == 0) then
295            write(*,*) "calchim: Error; no O1D tracer !!!"
[2302]296            call abort_physic("calchim","missing o1d tracer",1)
[1495]297         else
298            nbq = nbq + 1
299            niq(nbq) = i_o1d
300         end if
301         i_o2 = igcm_o2
302         if (i_o2 == 0) then
303            write(*,*) "calchim: Error; no O2 tracer !!!"
[2302]304            call abort_physic("calchim","missing o2 tracer",1)
[1495]305         else
306            nbq = nbq + 1
307            niq(nbq) = i_o2
308         end if
309         i_o3 = igcm_o3
310         if (i_o3 == 0) then
311            write(*,*) "calchim: Error; no O3 tracer !!!"
[2302]312            call abort_physic("calchim","missing o3 tracer",1)
[1495]313         else
314            nbq = nbq + 1
315            niq(nbq) = i_o3
316         end if
317         i_h = igcm_h
318         if (i_h == 0) then
319            write(*,*) "calchim: Error; no H tracer !!!"
[2302]320            call abort_physic("calchim","missing h tracer",1)
[1495]321         else
322            nbq = nbq + 1
323            niq(nbq) = i_h
324         end if
325         i_h2 = igcm_h2
326         if (i_h2 == 0) then
327            write(*,*) "calchim: Error; no H2 tracer !!!"
[2302]328            call abort_physic("calchim","missing h2 tracer",1)
[1495]329         else
330            nbq = nbq + 1
331            niq(nbq) = i_h2
332         end if
333         i_oh = igcm_oh
334         if (i_oh == 0) then
335            write(*,*) "calchim: Error; no OH tracer !!!"
[2302]336            call abort_physic("calchim","missing oh tracer",1)
[1495]337         else
338            nbq = nbq + 1
339            niq(nbq) = i_oh
340         end if
341         i_ho2 = igcm_ho2
342         if (i_ho2 == 0) then
343            write(*,*) "calchim: Error; no HO2 tracer !!!"
[2302]344            call abort_physic("calchim","missing ho2 tracer",1)
[1495]345         else
346            nbq = nbq + 1
347            niq(nbq) = i_ho2
348         end if
349         i_h2o2 = igcm_h2o2
350         if (i_h2o2 == 0) then
351            write(*,*) "calchim: Error; no H2O2 tracer !!!"
[2302]352            call abort_physic("calchim","missing h2o2 tracer",1)
[1495]353         else
354            nbq = nbq + 1
355            niq(nbq) = i_h2o2
356         end if
357         i_ch4 = igcm_ch4
358         if (i_ch4 == 0) then
359            write(*,*) "calchim: Error; no CH4 tracer !!!"
360            write(*,*) "CH4 will be ignored in the chemistry"
361         else
362            nbq = nbq + 1
363            niq(nbq) = i_ch4
364         end if
365         i_n2 = igcm_n2
366         if (i_n2 == 0) then
367            write(*,*) "calchim: Error; no N2 tracer !!!"
[2302]368            call abort_physic("calchim","missing n2 tracer",1)
[1495]369         else
370            nbq = nbq + 1
371            niq(nbq) = i_n2
372         end if
373         i_n = igcm_n
374         if (i_n == 0) then
375            if (photochem) then
376               write(*,*) "calchim: Error; no N tracer !!!"
[2302]377               call abort_physic("calchim","missing n tracer",1)
[1495]378            end if
379         else
380            nbq = nbq + 1
381            niq(nbq) = i_n
382         end if
383         i_n2d = igcm_n2d
384         if (i_n2d == 0) then
385            if (photochem) then
386               write(*,*) "calchim: Error; no N2D tracer !!!"
[2302]387               call abort_physic("calchim","missing n2d tracer",1)
[1495]388            end if
389         else
390            nbq = nbq + 1
391            niq(nbq) = i_n2d
392         end if
393         i_no = igcm_no
394         if (i_no == 0) then
395            if (photochem) then
396               write(*,*) "calchim: Error; no NO tracer !!!"
[2302]397               call abort_physic("calchim","missing no tracer",1)
[1495]398            end if
399         else
400            nbq = nbq + 1
401            niq(nbq) = i_no
402         end if
403         i_no2 = igcm_no2
404         if (i_no2 == 0) then
405            if (photochem) then
406               write(*,*) "calchim: Error; no NO2 tracer !!!"
[2302]407               call abort_physic("calchim","missing no2 tracer",1)
[1495]408            end if
409         else
410            nbq = nbq + 1
411            niq(nbq) = i_no2
412         end if
413         i_h2o = igcm_h2o_vap
414         if (i_h2o == 0) then
415            write(*,*) "calchim: Error; no water vapor tracer !!!"
[2302]416            call abort_physic("calchim","missing h2o_vap tracer",1)
[1495]417         else
418            nbq = nbq + 1
419            niq(nbq) = i_h2o
420         end if
[2461]421         i_hdo=igcm_hdo_vap
422         if (i_hdo == 0) then
423            write(*,*) "calchim: no HDO tracer !!!"
424!            call abort_physic("calchim","missing hdo_vap tracer",1)
425            write(*,*) "No deuterium chemistry considered"
426         else
427            nbq = nbq + 1
428            niq(nbq) = i_hdo
429            deutchem = .true.
430            write(*,*) "calchim: HDO tracer found in traceur.def"
431            write(*,*) "Deuterium chemistry included"
432         end if
433         i_od=igcm_od
434         if(deutchem) then
435            if (i_od == 0) then
436               write(*,*) "calchim: Error, no OD tracer !!!"
437               write(*,*) "OD is needed if HDO is in traceur.def"
438               call abort_physic("calchim","missing od tracer",1)
439            else
440               nbq = nbq + 1
441               niq(nbq) = i_od
442            end if
443         else
444            if (i_oplus /= 0) then
445               write(*,*) "calchim: Error: OD is present, but HDO is not!!!"
446               write(*,*) "Both must be in traceur.def if deuterium chemistry wanted"
447               call abort_physic("calchim","missing hdo tracer",1)
448            endif
449         endif
450         i_d=igcm_d
451         if(deutchem) then
452            if (i_d == 0) then
453               write(*,*) "calchim: Error, no D tracer !!!"
454               write(*,*) "D is needed if HDO is in traceur.def"
455               call abort_physic("calchim","missing d tracer",1)
456            else
457               nbq = nbq + 1
458               niq(nbq) = i_d
459            end if
460         else
461            if (i_d /= 0) then
462               write(*,*) "calchim: Error: D is present, but HDO is not!!!"
463               write(*,*) "Both must be in traceur.def if deuterium chemistry wanted"
464               call abort_physic("calchim","missing hdo tracer",1)
465            endif
466         endif
467         i_hd=igcm_hd
468         if(deutchem) then
469            if (i_hd == 0) then
470               write(*,*) "calchim: Error, no HD tracer !!!"
471               write(*,*) "HD is needed if HDO is in traceur.def"
472               call abort_physic("calchim","missing hd tracer",1)
473            else
474               nbq = nbq + 1
475               niq(nbq) = i_hd
476            end if
477         else
478            if (i_hd /= 0) then
479               write(*,*) "calchim: Error: HD is present, but HDO is not!!!"
480               write(*,*) "Both must be in traceur.def if deuterium chemistry wanted"
481               call abort_physic("calchim","missing hd tracer",1)
482            endif
483         endif
484         i_do2=igcm_do2
485         if(deutchem) then
486            if (i_do2 == 0) then
487               write(*,*) "calchim: Error, no DO2 tracer !!!"
488               write(*,*) "DO2 is needed if HDO is in traceur.def"
489               call abort_physic("calchim","missing do2 tracer",1)
490            else
491               nbq = nbq + 1
492               niq(nbq) = i_do2
493            end if
494         else
495            if (i_do2 /= 0) then
496               write(*,*) "calchim: Error: DO2 is present, but HDO is not!!!"
497               write(*,*) "Both must be in traceur.def if deuterium chemistry wanted"
498               call abort_physic("calchim","missing do2 tracer",1)
499            endif
500         endif
501         i_hdo2=igcm_hdo2
502         if(deutchem) then
503            if (i_hdo2 == 0) then
504               write(*,*) "calchim: Error, no HDO2 tracer !!!"
505               write(*,*) "HDO2 is needed if HDO is in traceur.def"
506               call abort_physic("calchim","missing hdo2 tracer",1)
507            else
508               nbq = nbq + 1
509               niq(nbq) = i_hdo2
510            end if
511         else
512            if (i_hdo2 /= 0) then
513               write(*,*) "calchim: Error: HDO2 is present, but HDO is not!!!"
514               write(*,*) "Both must be in traceur.def if deuterium chemistry wanted"
515               call abort_physic("calchim","missing hdo2 tracer",1)
516            endif
517         endif
[2158]518         i_o2plus = igcm_o2plus
519         if (i_o2plus == 0) then
[2170]520            write(*,*) "calchim: no O2+ tracer "
521            write(*,*) "Only neutral chemistry"
[2158]522         else
523            nbq = nbq + 1
524            niq(nbq) = i_o2plus
[2170]525            ionchem = .true.
526            write(*,*) "calchim: O2+ tracer found in traceur.def"
527            write(*,*) "Ion chemistry included"
[2158]528         end if
529         i_co2plus = igcm_co2plus
[2170]530         if(ionchem) then
531            if (i_co2plus == 0) then
532               write(*,*) "calchim: Error, no CO2+ tracer !!!"
533               write(*,*) "CO2+ is needed if O2+ is in traceur.def"
[2302]534               call abort_physic("calchim","missing co2plus tracer",1)
[2170]535            else
536               nbq = nbq + 1
537               niq(nbq) = i_co2plus
538            end if
[2158]539         else
[2170]540            if (i_co2plus /= 0) then
541               write(*,*) "calchim: Error: CO2+ is present, but O2+ is not!!!"
542               write(*,*) "Both must be in traceur.def if ion chemistry wanted"
[2302]543               call abort_physic("calchim","missing o2plus tracer",1)
[2170]544            endif
545         endif
[2158]546         i_oplus=igcm_oplus
[2170]547         if(ionchem) then
548            if (i_oplus == 0) then
549               write(*,*) "calchim: Error, no O+ tracer !!!"
550               write(*,*) "O+ is needed if O2+ is in traceur.def"
[2302]551               call abort_physic("calchim","missing oplus tracer",1)
[2170]552            else
553               nbq = nbq + 1
554               niq(nbq) = i_oplus
555            end if
[2158]556         else
[2170]557            if (i_oplus /= 0) then
558               write(*,*) "calchim: Error: O+ is present, but O2+ is not!!!"
559               write(*,*) "Both must be in traceur.def if ion chemistry wanted"
[2302]560               call abort_physic("calchim","missing o2plus tracer",1)
[2170]561            endif
562         endif
[2158]563         i_noplus=igcm_noplus
[2170]564         if(ionchem) then
565            if (i_noplus == 0) then
566               write(*,*) "calchim: Error, no NO+ tracer !!!"
567               write(*,*) "NO+ is needed if O2+ is in traceur.def"
[2302]568               call abort_physic("calchim","missing noplus tracer",1)
[2170]569            else
570               nbq = nbq + 1
571               niq(nbq) = i_noplus
572            end if
[2158]573         else
[2170]574            if (i_noplus /= 0) then
575               write(*,*) "calchim: Error: NO+ is present, but O2+ is not!!!"
576               write(*,*) "Both must be in traceur.def if ion chemistry wanted"
577            endif
578         endif
[2158]579         i_coplus=igcm_coplus
[2170]580         if(ionchem) then
581            if (i_coplus == 0) then
582               write(*,*) "calchim: Error, no CO+ tracer !!!"
583               write(*,*) "CO+ is needed if O2+ is in traceur.def"
[2302]584               call abort_physic("calchim","missing coplus tracer",1)               
[2170]585            else
586               nbq = nbq + 1
587               niq(nbq) = i_coplus
588            end if
[2158]589         else
[2170]590            if (i_coplus /= 0) then
591               write(*,*) "calchim: Error: CO+ is present, but O2+ is not!!!"
592               write(*,*) "Both must be in traceur.def if ion chemistry wanted"
593            endif
594         endif
[2158]595         i_cplus=igcm_cplus
[2170]596         if(ionchem) then
597            if (i_cplus == 0) then
598               write(*,*) "calchim: Error, no C+ tracer !!!"
599               write(*,*) "C+ is needed if O2+ is in traceur.def"
[2302]600               call abort_physic("calchim","missing cplus tracer",1)
[2170]601            else
602               nbq = nbq + 1
603               niq(nbq) = i_cplus
604            end if
[2158]605         else
[2170]606            if (i_cplus /= 0) then
607               write(*,*) "calchim: Error: C+ is present, but O2+ is not!!!"
608               write(*,*) "Both must be in traceur.def if ion chemistry wanted"
609            endif
610         endif
[2158]611         i_n2plus=igcm_n2plus
[2170]612         if(ionchem) then
613            if (i_n2plus == 0) then
614               write(*,*) "calchim: Error, no N2+ tracer !!!"
615               write(*,*) "N2+ is needed if O2+ is in traceur.def"
[2302]616               call abort_physic("calchim","missing n2plus tracer",1)
[2170]617            else
618               nbq = nbq + 1
619               niq(nbq) = i_n2plus
620            end if
[2158]621         else
[2170]622            if (i_n2plus /= 0) then
623               write(*,*) "calchim: Error: N2+ is present, but O2+ is not!!!"
624               write(*,*) "Both must be in traceur.def if ion chemistry wanted"
625            endif
626         endif
[2158]627         i_nplus=igcm_nplus
[2170]628         if(ionchem) then
629            if (i_nplus == 0) then
630               write(*,*) "calchim: Error, no N+ tracer !!!"
631               write(*,*) "N+ is needed if O2+ is in traceur.def"
[2302]632               call abort_physic("calchim","missing nplus tracer",1)
[2170]633            else
634               nbq = nbq + 1
635               niq(nbq) = i_nplus
636            end if
[2158]637         else
[2170]638            if (i_nplus /= 0) then
639               write(*,*) "calchim: Error: N+ is present, but O2+ is not!!!"
640               write(*,*) "Both must be in traceur.def if ion chemistry wanted"
641            endif
642         endif
[2158]643         i_hplus=igcm_hplus
[2170]644         if(ionchem) then
645            if (i_hplus == 0) then
646               write(*,*) "calchim: Error, no H+ tracer !!!"
647               write(*,*) "H+ is needed if O2+ is in traceur.def"
[2302]648               call abort_physic("calchim","missing hplus tracer",1)
[2170]649            else
650               nbq = nbq + 1
651               niq(nbq) = i_hplus
652            end if
[2158]653         else
[2170]654            if (i_hplus /= 0) then
655               write(*,*) "calchim: Error: H+ is present, but O2+ is not!!!"
656               write(*,*) "Both must be in traceur.def if ion chemistry wanted"
657            endif
658         endif
[2158]659         i_hco2plus=igcm_hco2plus
[2170]660         if(ionchem) then
661            if (i_hco2plus == 0) then
662               write(*,*) "calchim: Error, no HCO2+ tracer !!!"
663               write(*,*) "HCO2+ is needed if O2+ is in traceur.def"
[2302]664               call abort_physic("calchim","missing hco2plus tracer",1)
[2170]665            else
666               nbq = nbq + 1
667               niq(nbq) = i_hco2plus
668            end if
[2158]669         else
[2170]670            if (i_hco2plus /= 0) then
671               write(*,*) "calchim: Error: HCO2+ is present, but O2+ is not!!!"
672               write(*,*) "Both must be in traceur.def if ion chemistry wanted"
673            endif
674         endif
[2284]675         i_hcoplus=igcm_hcoplus
676         if(ionchem) then
677            if (i_hcoplus == 0) then
678               write(*,*) "calchim: Error, no HCO+ tracer !!!"
679               write(*,*) "HCO+ is needed if O2+ is in traceur.def"
[2302]680               call abort_physic("calchim","missing hcoplus tracer",1)
[2284]681            else
682               nbq = nbq + 1
683               niq(nbq) = i_hcoplus
684            end if
685         else
686            if (i_hcoplus /= 0) then
687               write(*,*) "calchim: Error: HCO+ is present, but O2+ is not!!!"
688               write(*,*) "Both must be in traceur.def if ion chemistry wanted"
689            endif
690         endif
[2302]691         i_h2oplus=igcm_h2oplus
692         if(ionchem) then
693            if (i_h2oplus == 0) then
694               write(*,*) "calchim: Error, no H2O+ tracer !!!"
695               write(*,*) "H2O+ is needed if O2+ is in traceur.def"
696               call abort_physic("calchim","missing h2oplus tracer",1)
697            else
698               nbq = nbq + 1
699               niq(nbq) = i_h2oplus
700            end if
701         else
702            if (i_h2oplus /= 0) then
703               write(*,*) "calchim: Error: H2O+ is present, but O2+ is not!!!"
704               write(*,*) "Both must be in traceur.def if ion chemistry wanted"
705            endif
706         endif
[2321]707         i_h3oplus=igcm_h3oplus
708         if(ionchem) then
709            if (i_h3oplus == 0) then
710               write(*,*) "calchim: Error, no H3O+ tracer !!!"
711               write(*,*) "H3O+ is needed if O2+ is in traceur.def"
712               call abort_physic("calchim","missing h3oplus tracer",1)
713            else
714               nbq = nbq + 1
715               niq(nbq) = i_h3oplus
716            end if
717         else
718            if (i_h3oplus /= 0) then
719               write(*,*) "calchim: Error: H3O+ is present, but O2+ is not!!!"
720               write(*,*) "Both must be in traceur.def if ion chemistry wanted"
721            endif
722         endif
723         i_ohplus=igcm_ohplus
724         if(ionchem) then
725            if (i_ohplus == 0) then
726               write(*,*) "calchim: Error, no OH+ tracer !!!"
727               write(*,*) "OH+ is needed if O2+ is in traceur.def"
728               call abort_physic("calchim","missing ohplus tracer",1)
729            else
730               nbq = nbq + 1
731               niq(nbq) = i_ohplus
732            end if
733         else
734            if (i_ohplus /= 0) then
735               write(*,*) "calchim: Error: OH+ is present, but O2+ is not!!!"
736               write(*,*) "Both must be in traceur.def if ion chemistry wanted"
737            endif
738         endif
[2158]739         i_elec = igcm_elec
[2170]740         if(ionchem) then
741            if (i_elec == 0) then
742               write(*,*) "calchim: Error, no e- tracer !!!"
743               write(*,*) "e- is needed if O2+ is in traceur.def"
[2302]744               call abort_physic("calchim","missing elec tracer",1)
[2170]745            else
746               nbq = nbq + 1
747               niq(nbq) = i_elec
748            end if
[2158]749         else
[2170]750            if (i_elec /= 0) then
751               write(*,*) "calchim: Error: e- is present, but O2+ is not!!!"
752               write(*,*) "Both must be in traceur.def if ion chemistry wanted"
753            endif
754         endif
[1495]755         write(*,*) 'calchim: found nbq    = ',nbq,' tracers'
[2461]756
757         write(*,*) 'calchim: tracer indices=',niq(:)
758
759
760         if (photochem) then
761            if (jonline) then
762               print*,'calchim: Read UV absorption cross-sections'
763               !Add two photodissociations in deuterium chemistry included
764               if(deutchem) nphot = nphot + 2
765               call init_photolysis     ! on-line photolysis
766            else
767               print*,'calchim: Read photolysis lookup table'
768               call read_phototable     ! off-line photolysis
769            end if
770         end if
771
772         if(.not.unichim) then
773            !Read reaction rates from external file if the upper atmospheric
774            !chemistry is called
775            call chemthermos_readini
776         endif
777
778!         stop
[1495]779         firstcall = .false.
780      end if ! if (firstcall)
781
782! Initializations
783
784      zycol(:,:)    = 0.
785      dqchim(:,:,:) = 0.
786      dqschim(:,:)  = 0.
787
[2030]788      kb = 1.3806e-23
789
[1495]790!     latvl1= 22.27
791!     lonvl1= -47.94
[2029]792!     ig_vl1= 1+ int( (1.5-(latvl1-90.)*48./180.)  -2 )*64. +    &
793!             int(1.5+(lonvl1+180)*64./360.)
[1495]794
795!=======================================================================
796!     loop over grid
797!=======================================================================
798
799      do ig = 1,ngrid
800         
801         foundswitch = 0
802         do l = 1,nlayer
803            do i = 1,nbq
804               iq = niq(i) ! get tracer index
805               zq(ig,l,iq) = pq(ig,l,iq) + pdq(ig,l,iq)*ptimestep
806               zycol(l,iq) = zq(ig,l,iq)*mmean(ig,l)/mmol(iq)
807            end do
808            zt(ig,l)  = pt(ig,l) + pdt(ig,l)*ptimestep
809            zu(ig,l)  = pu(ig,l) + pdu(ig,l)*ptimestep
810            zv(ig,l)  = pv(ig,l) + pdv(ig,l)*ptimestep
811            zpress(l) = pplay(ig,l)/100.
812            ztemp(l)  = zt(ig,l)
813            zdens(l)  = zpress(l)/(kb*1.e4*ztemp(l))
814            zlocal(l) = zzlay(ig,l)/1000.
815            zmmean(l) = mmean(ig,l)
[2158]816            ztelec(l) = temp_elect(zlocal(l),ztemp(l),1)
817            !Electronic temperature. Index 1 -> Viking; Index 2-> MAVEN
[1495]818
819!           surfdust1d and surfice1d: conversion from m2/m3 to cm2/cm3
820
821            surfdust1d(l) = surfdust(ig,l)*1.e-2
822            surfice1d(l)  = surfice(ig,l)*1.e-2
823
[2158]824!           search for switch index between regions 
[2170]825
826            if (unichim) then
827               lswitch = nlayer + 1
[2158]828            else
[2461]829               if (foundswitch == 0 .and. pplay(ig,l) < 10.) then
[1495]830                  lswitch = l
831                  foundswitch = 1
832               end if
[2158]833            endif
[1495]834         end do ! of do l=1,nlayer
835
836         szacol = acos(mu0(ig))*180./pi
[2031]837         taucol = tau(ig)
[1495]838
839!=======================================================================
840!     call chemical subroutines
841!=======================================================================
842
843         if (photochem) then
[2170]844            ! set number of reactions, depending on ion chemistry or not
[2461]845            nb_reaction_4_ion  = 64
846            nb_reaction_4_deut = 35
847
848            !Default numbers if no ion and no deuterium chemistry included
849
850            nb_reaction_4_max = 31     ! set number of bimolecular reactions
851            nb_reaction_3_max = 6      ! set number of quadratic reactions
852            nquench           = 9      ! set number of quenching + heterogeneous
853            nphotion          = 0      ! set number of photoionizations
854
[2170]855            if (ionchem) then
[2461]856               nb_reaction_4_max = nb_reaction_4_max + nb_reaction_4_ion 
[2170]857               nphotion          = 18   ! set number of photoionizations
[2461]858            endif
859            if(deutchem) then
860               nb_reaction_4_max = nb_reaction_4_max + nb_reaction_4_deut 
[2170]861            end if
862
863!        nb_phot_max is the total number of processes that are treated
864!        numerically as a photolysis:
865
866            nb_phot_max = nphot + nphotion + nquench
867
868!        call main photochemistry routine
869
[2461]870            call photochemistry(nlayer,nq,nbq,ionchem,deutchem,           &
871                           nb_reaction_3_max,nb_reaction_4_max,nphot,     &
872                           nb_phot_max,nphotion,                          &
[2433]873                           jonline,ig,lswitch,zycol,szacol,ptimestep,     &
[2170]874                           zpress,zlocal,ztemp,ztelec,zdens,zmmean,       &
875                           dist_sol,zday,surfdust1d,surfice1d,            &
[2321]876                           jo3,jh2o,em_no,em_o2,taucol,iter)
[1495]877
878!        ozone photolysis, for output
879
880            do l = 1,nlayer
881               jo3_3d(ig,l) = jo3(l)
[2030]882               jh2o_3d(ig,l) = jh2o(l)
[1495]883               iter_3d(ig,l) = iter(l)
884            end do
[2170]885
[1495]886!        condensation of h2o2
887
888            call perosat(ngrid, nlayer, nq,                        &
889                         ig,ptimestep,pplev,pplay,                 &
890                         ztemp,zycol,dqcloud,dqscloud)
891
[2170]892!        case of separate photochemical model in the thermosphere
893
894            if (.not.unichim) then
895               chemthermod = 3   !C/O/H/N/ions
[2158]896               call chemthermos(ig,nlayer,lswitch,chemthermod,zycol,ztemp,&
897                    zdens,zpress,zlocal,szacol,ptimestep,zday,&
898                    em_no,em_o2)
899            end if
[1495]900
[2321]901            do l = 1,nlayer
902               emission_no(ig,l) = em_no(l)
903               emission_o2(ig,l) = em_o2(l)
904            end do
[2170]905         end if  ! photochem
906
[1495]907!        dry deposition
908
909         if (depos) then
910            call deposition(ngrid, nlayer, nq,                      &
911                            ig, ig_vl1, pplay, pplev, zzlay, zzlev, &
912                            zu, zv, zt, zycol, ptimestep, co2ice)
913         end if
914
915!=======================================================================
916!     tendencies
917!=======================================================================
918
919!     index of the most abundant species at each level
920
921!         major(:) = maxloc(zycol, dim = 2)
922
923!     tendency for the most abundant species = - sum of others
924
925         do l = 1,nlayer
926            iloc=maxloc(zycol(l,:))
927            iqmax=iloc(1)
928            do i = 1,nbq
929               iq = niq(i) ! get tracer index
930               if (iq /= iqmax) then
931                  dqchim(ig,l,iq) = (zycol(l,iq)*mmol(iq)/mmean(ig,l)  &
932                                   - zq(ig,l,iq))/ptimestep
933                  dqchim(ig,l,iqmax) = dqchim(ig,l,iqmax)              &
934                                     - dqchim(ig,l,iq)
935               end if
936            end do
937         end do ! of do l = 1,nlayer
938
939!=======================================================================
940!     end of loop over grid
941!=======================================================================
942
943      end do ! of do ig=1,ngrid
944
945!=======================================================================
946!     write outputs
947!=======================================================================
948
949! value of parameter 'output' to trigger writting of outputs
950! is set above at the declaration of the variable.
951
952      if (photochem .and. output) then
953         if (ngrid > 1) then
954            call writediagfi(ngrid,'jo3','j o3->o1d',    &
955                             's-1',3,jo3_3d(1,1))
[2030]956            call writediagfi(ngrid,'jh2o','jh2o',    &
957                             's-1',3,jh2o_3d(1,1))
[1495]958            call writediagfi(ngrid,'iter','iterations',  &
959                             ' ',3,iter_3d(1,1))
[2170]960
[2321]961!            if (.not. unichim) then
[2158]962               call writediagfi(ngrid,'emission_no',        &
963                    'NO nightglow emission rate','cm-3 s-1',3,emission_no)
964               call writediagfi(ngrid,'emission_o2',        &
965                    'O2 nightglow emission rate','cm-3 s-1',3,emission_o2)
[2321]966!            endif
[2042]967           
[1495]968               call wstats(ngrid,'jo3','j o3->o1d',       &
969                           's-1',3,jo3_3d(1,1))
[1888]970               call wstats(ngrid,'emission_no',           &
971                   'NO nightglow emission rate','cm-3 s-1',3,emission_no)
972               call wstats(ngrid,'emission_o2',           &
973                   'O2 nightglow emission rate','cm-3 s-1',3,emission_o2)
[1495]974               call wstats(ngrid,'mmean','mean molecular mass',       &
975                           'g.mole-1',3,mmean(1,1))
976         end if ! of if (ngrid.gt.1)
977      end if ! of if (output)
978
[2162]979      end subroutine calchim
980
981
982    subroutine ini_calchim_mod(ngrid,nlayer,nq)
983 
984      implicit none
985 
986      integer,intent(in) :: ngrid ! number of atmospheric columns
987      integer,intent(in) :: nlayer ! number of atmospheric layers
988      integer,intent(in) :: nq ! number of tracers
989
990      allocate(zdqchim(ngrid,nlayer,nq))
991      zdqchim(:,:,:)=0
992      allocate(zdqschim(ngrid,nq))
993      zdqschim(:,:)=0
994
995    end subroutine ini_calchim_mod
996
997
998    subroutine end_calchim_mod
999
1000      implicit none
1001
1002      if (allocated(zdqchim))      deallocate(zdqchim)
1003      if (allocated(zdqschim))      deallocate(zdqschim)
1004
1005    end subroutine end_calchim_mod
1006
1007END MODULE calchim_mod
1008
Note: See TracBrowser for help on using the repository browser.