source: trunk/LMDZ.MARS/libf/aeronomars/calchim.F90 @ 1198

Last change on this file since 1198 was 1125, checked in by flefevre, 11 years ago

calchim.F90 et photochemistry.F:

Suppression des references a nlayermx et nqmx.
=> Passage de nlayer et nq en arguments de subroutines.

File size: 26.6 KB
Line 
1      subroutine calchim(ngrid,nlayer,nq,                                   &
2                         ptimestep,pplay,pplev,pt,pdt,dist_sol,mu0,         &
3                         zzlev,zzlay,zday,pq,pdq,dqchim,dqschim,dqcloud,    &
4                         dqscloud,tauref,co2ice,                            &
5                         pu,pdu,pv,pdv,surfdust,surfice)
6
7      use tracer_mod, only: igcm_co2, igcm_co, igcm_o, igcm_o1d, igcm_o2, &
8                            igcm_o3, igcm_h, igcm_h2, igcm_oh, igcm_ho2,  &
9                            igcm_h2o2, igcm_ch4, igcm_n2, igcm_h2o_vap,   &
10                            igcm_no, igcm_n, igcm_no2, igcm_n2d,          &
11                            igcm_o2plus, igcm_co2plus, igcm_oplus,        &
12                            igcm_coplus, igcm_cplus, igcm_nplus,          &
13                            igcm_noplus, igcm_n2plus, igcm_hplus,         &
14                            igcm_hco2plus, igcm_elec, mmol
15
16      use conc_mod, only: mmean ! mean molecular mass of the atmosphere
17
18      implicit none
19
20!=======================================================================
21!
22!   subject:
23!   --------
24!
25!  Prepare the call for the photochemical module, and send back the
26!  tendencies from photochemistry in the chemical species mass mixing ratios
27!
28!   Author:   Sebastien Lebonnois (08/11/2002)
29!   -------
30!    update 12/06/2003 for water ice clouds and compatibility with dust
31!    update 07/2003 for coupling with thermosphere (Monica Angelats-i-Coll)
32!    update 03/05/2005 cosmetic changes (Franck Lefevre)
33!    update sept. 2008 identify tracers by their names (Ehouarn Millour)
34!    update 05/12/2011 synchronize with latest version of chemistry (Franck Lefevre)
35!    update 16/03/2012 optimization (Franck Lefevre)
36!    update 11/12/2013 optimization (Franck Lefevre)
37!
38!   Arguments:
39!   ----------
40!
41!  Input:
42!
43!    ptimestep                  timestep (s)
44!    pplay(ngrid,nlayer)    Pressure at the middle of the layers (Pa)
45!    pplev(ngrid,nlayer+1)  Intermediate pressure levels (Pa)
46!    pt(ngrid,nlayer)       Temperature (K)
47!    pdt(ngrid,nlayer)      Temperature tendency (K)
48!    pu(ngrid,nlayer)       u component of the wind (ms-1)
49!    pdu(ngrid,nlayer)      u component tendency (K)
50!    pv(ngrid,nlayer)       v component of the wind (ms-1)
51!    pdv(ngrid,nlayer)      v component tendency (K)
52!    dist_sol                   distance of the sun (AU)
53!    mu0(ngrid)               cos of solar zenith angle (=1 when sun at zenith)
54!    pq(ngrid,nlayer,nq)  Advected fields, ie chemical species here
55!    pdq(ngrid,nlayer,nq) Previous tendencies on pq
56!    tauref(ngrid)            Optical depth at 7 hPa
57!    co2ice(ngrid)            co2 ice surface layer (kg.m-2)
58!    surfdust(ngrid,nlayer) dust surface area (m2/m3)
59!    surfice(ngrid,nlayer)  ice surface area (m2/m3)
60!
61!  Output:
62!
63!    dqchim(ngrid,nlayer,nq) ! tendencies on pq due to chemistry
64!    dqschim(ngrid,nq)         ! tendencies on qsurf
65!
66!=======================================================================
67
68!#include "dimensions.h"
69!#include "dimphys.h"
70#include "chimiedata.h"
71!#include "tracer.h"
72#include "comcstfi.h"
73#include "callkeys.h"
74!#include "conc.h"
75
76!     input:
77
78      integer,intent(in) :: ngrid    ! number of atmospheric columns
79      integer,intent(in) :: nlayer   ! number of atmospheric layers
80      integer,intent(in) :: nq       ! number of tracers
81      real :: ptimestep
82      real :: pplay(ngrid,nlayer)    ! pressure at the middle of the layers
83      real :: zzlay(ngrid,nlayer)    ! pressure at the middle of the layers
84      real :: pplev(ngrid,nlayer+1)  ! intermediate pressure levels
85      real :: zzlev(ngrid,nlayer+1)  ! altitude at layer boundaries
86      real :: pt(ngrid,nlayer)       ! temperature
87      real :: pdt(ngrid,nlayer)      ! temperature tendency
88      real :: pu(ngrid,nlayer)       ! u component of the wind (m.s-1)
89      real :: pdu(ngrid,nlayer)      ! u component tendency
90      real :: pv(ngrid,nlayer)       ! v component of the wind (m.s-1)
91      real :: pdv(ngrid,nlayer)      ! v component tendency
92      real :: dist_sol               ! distance of the sun (AU)
93      real :: mu0(ngrid)             ! cos of solar zenith angle (=1 when sun at zenith)
94      real :: pq(ngrid,nlayer,nq)    ! tracers mass mixing ratio
95      real :: pdq(ngrid,nlayer,nq)   ! previous tendencies
96      real :: zday                   ! date (time since Ls=0, in martian days)
97      real :: tauref(ngrid)          ! optical depth at 7 hPa
98      real :: co2ice(ngrid)          ! co2 ice surface layer (kg.m-2)
99      real :: surfdust(ngrid,nlayer) ! dust surface area (m2/m3)
100      real :: surfice(ngrid,nlayer)  !  ice surface area (m2/m3)
101
102!     output:
103
104      real :: dqchim(ngrid,nlayer,nq)   ! tendencies on pq due to chemistry
105      real :: dqschim(ngrid,nq)         ! tendencies on qsurf
106      real :: dqcloud(ngrid,nlayer,nq)  ! tendencies on pq due to condensation
107      real :: dqscloud(ngrid,nq)        ! tendencies on qsurf
108
109!     local variables:
110
111      integer,save :: nbq                   ! number of tracers used in the chemistry
112      integer,allocatable,save :: niq(:)    ! array storing the indexes of the tracers
113      integer :: iloc(1)                    ! index of major species
114      integer :: ig,l,i,iq,iqmax
115      integer :: foundswitch, lswitch
116      integer,save :: chemthermod
117
118      integer,save :: i_co2  = 0
119      integer,save :: i_co   = 0
120      integer,save :: i_o    = 0
121      integer,save :: i_o1d  = 0
122      integer,save :: i_o2   = 0
123      integer,save :: i_o3   = 0
124      integer,save :: i_h    = 0
125      integer,save :: i_h2   = 0
126      integer,save :: i_oh   = 0
127      integer,save :: i_ho2  = 0
128      integer,save :: i_h2o2 = 0
129      integer,save :: i_ch4  = 0
130      integer,save :: i_n2   = 0
131      integer,save :: i_h2o  = 0
132      integer,save :: i_n    = 0
133      integer,save :: i_no   = 0
134      integer,save :: i_no2  = 0
135      integer,save :: i_n2d  = 0
136      integer,save :: i_co2plus=0
137      integer,save :: i_oplus=0
138      integer,save :: i_o2plus=0
139      integer,save :: i_coplus=0
140      integer,save :: i_cplus=0
141      integer,save :: i_nplus=0
142      integer,save :: i_noplus=0
143      integer,save :: i_n2plus=0
144      integer,save :: i_hplus=0
145      integer,save :: i_hco2plus=0
146      integer,save :: i_elec=0
147
148      integer :: ig_vl1
149
150      real    :: latvl1, lonvl1
151      real    :: zq(ngrid,nlayer,nq)   ! pq+pdq*ptimestep before chemistry
152                                       ! new mole fraction after
153      real    :: zt(ngrid,nlayer)      ! temperature
154      real    :: zu(ngrid,nlayer)      ! u component of the wind
155      real    :: zv(ngrid,nlayer)      ! v component of the wind
156      real    :: taucol                ! optical depth at 7 hPa
157
158      logical,save :: firstcall = .true.
159      logical,save :: depos = .false.  ! switch for dry deposition
160
161!     for each column of atmosphere:
162
163      real :: zpress(nlayer)       ! Pressure (mbar)
164      real :: zdens(nlayer)        ! Density  (cm-3)
165      real :: ztemp(nlayer)        ! Temperature (K)
166      real :: zlocal(nlayer)       ! Altitude (km)
167      real :: zycol(nlayer,nq)     ! Composition (mole fractions)
168      real :: szacol               ! Solar zenith angle
169      real :: surfice1d(nlayer)    ! Ice surface area (cm2/cm3)
170      real :: surfdust1d(nlayer)   ! Dust surface area (cm2/cm3)
171      real :: jo3(nlayer)          ! Photodissociation rate O3->O1D (s-1)
172
173!     for output:
174
175      logical :: output             ! to issue calls to writediagfi and stats
176      parameter (output = .true.)
177      real :: jo3_3d(ngrid,nlayer)  ! Photodissociation rate O3->O1D (s-1)
178
179!=======================================================================
180!     initialization of the chemistry (first call only)
181!=======================================================================
182
183      if (firstcall) then
184
185         if (photochem) then
186            print*,'calchim: Read photolysis lookup table'
187            call read_phototable
188         end if
189         ! find index of chemical tracers to use
190         allocate(niq(nq))
191         ! Listed here are all tracers that can go into photochemistry
192         nbq = 0 ! to count number of tracers
193         ! Species ALWAYS present if photochem=.T. or thermochem=.T.
194         i_co2 = igcm_co2
195         if (i_co2 == 0) then
196            write(*,*) "calchim: Error; no CO2 tracer !!!"
197            stop
198         else
199            nbq = nbq + 1
200            niq(nbq) = i_co2
201         end if
202         i_co = igcm_co
203         if (i_co == 0) then
204            write(*,*) "calchim: Error; no CO tracer !!!"
205            stop
206         else
207            nbq = nbq + 1
208            niq(nbq) = i_co
209         end if
210         i_o = igcm_o
211         if (i_o == 0) then
212            write(*,*) "calchim: Error; no O tracer !!!"
213            stop
214         else
215            nbq = nbq + 1
216            niq(nbq) = i_o
217         end if
218         i_o1d = igcm_o1d
219         if (i_o1d == 0) then
220            write(*,*) "calchim: Error; no O1D tracer !!!"
221            stop
222         else
223            nbq = nbq + 1
224            niq(nbq) = i_o1d
225         end if
226         i_o2 = igcm_o2
227         if (i_o2 == 0) then
228            write(*,*) "calchim: Error; no O2 tracer !!!"
229            stop
230         else
231            nbq = nbq + 1
232            niq(nbq) = i_o2
233         end if
234         i_o3 = igcm_o3
235         if (i_o3 == 0) then
236            write(*,*) "calchim: Error; no O3 tracer !!!"
237            stop
238         else
239            nbq = nbq + 1
240            niq(nbq) = i_o3
241         end if
242         i_h = igcm_h
243         if (i_h == 0) then
244            write(*,*) "calchim: Error; no H tracer !!!"
245            stop
246         else
247            nbq = nbq + 1
248            niq(nbq) = i_h
249         end if
250         i_h2 = igcm_h2
251         if (i_h2 == 0) then
252            write(*,*) "calchim: Error; no H2 tracer !!!"
253            stop
254         else
255            nbq = nbq + 1
256            niq(nbq) = i_h2
257         end if
258         i_oh = igcm_oh
259         if (i_oh == 0) then
260            write(*,*) "calchim: Error; no OH tracer !!!"
261            stop
262         else
263            nbq = nbq + 1
264            niq(nbq) = i_oh
265         end if
266         i_ho2 = igcm_ho2
267         if (i_ho2 == 0) then
268            write(*,*) "calchim: Error; no HO2 tracer !!!"
269            stop
270         else
271            nbq = nbq + 1
272            niq(nbq) = i_ho2
273         end if
274         i_h2o2 = igcm_h2o2
275         if (i_h2o2 == 0) then
276            write(*,*) "calchim: Error; no H2O2 tracer !!!"
277            stop
278         else
279            nbq = nbq + 1
280            niq(nbq) = i_h2o2
281         end if
282         i_ch4 = igcm_ch4
283         if (i_ch4 == 0) then
284            write(*,*) "calchim: Error; no CH4 tracer !!!"
285            write(*,*) "CH4 will be ignored in the chemistry"
286         else
287            nbq = nbq + 1
288            niq(nbq) = i_ch4
289         end if
290         i_n2 = igcm_n2
291         if (i_n2 == 0) then
292            write(*,*) "calchim: Error; no N2 tracer !!!"
293            stop
294         else
295            nbq = nbq + 1
296            niq(nbq) = i_n2
297         end if
298         i_h2o = igcm_h2o_vap
299         if (i_h2o == 0) then
300            write(*,*) "calchim: Error; no water vapor tracer !!!"
301            stop
302         else
303            nbq = nbq + 1
304            niq(nbq) = i_h2o
305         end if
306         !Check tracers needed for thermospheric chemistry
307         if(thermochem) then
308            chemthermod=0  !Default: C/O/H chemistry
309            !Nitrogen chemistry
310            !NO is used to determine if N chemistry is wanted
311            !chemthermod=2 -> N chemistry
312            i_no = igcm_no
313            if (i_no == 0) then
314               write(*,*) "calchim: no NO tracer"
315               write(*,*) "C/O/H themosp chemistry only "
316            else
317               nbq = nbq + 1
318               niq(nbq) = i_no
319               chemthermod=2
320               write(*,*) "calchim: NO in traceur.def"
321               write(*,*) "Nitrogen chemistry included"
322            end if
323            ! N
324            i_n = igcm_n
325            if(chemthermod == 2) then
326               if (i_n == 0) then
327                  write(*,*) "calchim: Error; no N tracer !!!"
328                  write(*,*) "N is needed if NO is in traceur.def"
329                  stop
330               else
331                  nbq = nbq + 1
332                  niq(nbq) = i_n
333               end if
334            else
335               if (i_n /= 0) then
336                  write(*,*) "calchim: Error: N is present, but NO is not!!!"
337                  write(*,*) "Both must be in traceur.def if N chemistry is wanted"
338                  stop
339               endif
340            endif    !Of if(chemthermod == 2)
341            ! NO2
342            i_no2 = igcm_no2
343            if(chemthermod == 2) then
344               if (i_no2 == 0) then
345                  write(*,*) "calchim: Error; no NO2 tracer !!!"
346                  write(*,*) "NO2 is needed if NO is in traceur.def"
347                  stop
348               else
349                  nbq = nbq + 1
350                  niq(nbq) = i_no2
351               end if
352            else
353               if (i_no2 /= 0) then
354                  write(*,*) "calchim: Error: N is present, but NO is not!!!"
355                  write(*,*) "Both must be in traceur.def if N chemistry is wanted"
356                  stop
357               endif
358            endif     !Of if(chemthermod == 2)
359            ! N(2D)
360            if(chemthermod == 2) then
361               i_n2d = igcm_n2d
362               if (i_n2d == 0) then
363                  write(*,*) "calchim: Error; no N2D !!!"
364                  write(*,*) "N2D is needed if NO is in traceur.def"
365                  stop
366               else
367                  nbq = nbq + 1
368                  niq(nbq) = i_n2d
369               end if
370            else
371               if (i_n2d /= 0) then
372                  write(*,*) "calchim: Error: N2D is present, but NO is not!!!"
373                  write(*,*) "Both must be in traceur.def if N chemistry wanted"
374                  stop
375               endif
376            endif    !Of if(chemthermod == 2)
377            ! Ions
378            ! O2+ is used to determine if ion chemistry is needed
379            ! chemthermod=3 -> ion chemistry
380            i_o2plus = igcm_o2plus
381            if(chemthermod == 2) then
382               if (i_o2plus == 0) then
383                  write(*,*) "calchim: no O2+ tracer; no ion chemistry"
384               else
385                  nbq = nbq + 1
386                  niq(nbq) = i_o2plus
387                  chemthermod = 3
388                  write(*,*) "calchim: O2+ in traceur.def"
389                  write(*,*) "Ion chemistry included"
390               end if
391            else
392               if (i_o2plus /= 0) then
393                  write(*,*) "calchim: O2+ is present, but NO is not!!!"
394                  write(*,*) "Both must be in traceur.def if ionosphere wanted"
395                  stop
396               endif
397            endif   !Of if(chemthermod == 2)
398            ! CO2+
399            i_co2plus = igcm_co2plus
400            if(chemthermod == 3) then
401               if (i_co2plus == 0) then
402                  write(*,*) "calchim: Error; no CO2+ tracer !!!"
403                  write(*,*) "CO2+ is needed if O2+ is in traceur.def"
404                  stop
405               else
406                  nbq = nbq + 1
407                  niq(nbq) = i_co2plus
408               end if
409            else
410               if (i_co2plus /= 0) then
411                  write(*,*) "calchim: Error: CO2+ is present, but O2+ is not!!!"
412                  write(*,*) "Both must be in traceur.def if ionosphere wanted"
413                  stop
414               endif
415            endif    !Of if(chemthermod == 3)
416            ! O+
417            i_oplus = igcm_oplus
418            if(chemthermod == 3) then
419               if (i_oplus == 0) then
420                  write(*,*) "calchim: Error; no O+ tracer !!!"
421                  write(*,*) "O+ is needed if O2+ is in traceur.def"
422                  stop
423               else
424                  nbq = nbq + 1
425                  niq(nbq) = i_oplus
426               end if
427            else
428               if (i_oplus /= 0) then
429                  write(*,*) "calchim: Error: O+ is present, but O2+ is not!!!"
430                  write(*,*) "Both must be in traceur.def if ionosphere wanted"
431                  stop
432               endif
433            endif   !Of if (chemthermod == 3)
434            ! CO+
435            i_coplus = igcm_coplus
436            if(chemthermod == 3) then
437               if (i_coplus == 0) then
438                  write(*,*) "calchim: Error; no CO+ tracer !!!"
439                  write(*,*) "CO+ is needed if O2+ is in traceur.def"
440                  stop
441               else
442                  nbq = nbq + 1
443                  niq(nbq) = i_coplus
444               end if
445            else
446               if (i_coplus /= 0) then
447                  write(*,*) "calchim: Error: CO+ is present, but O2+ is not!!!"
448                  write(*,*) " Both must be in traceur.def if ionosphere wanted"
449                  stop
450               endif
451            endif   ! Of if (chemthermod == 3)
452            ! C+
453            i_cplus = igcm_cplus
454            if(chemthermod == 3) then
455               if (i_cplus == 0) then
456                  write(*,*) "calchim: Error; no C+ tracer !!!"
457                  write(*,*) "C+ is needed if O2+ is in traceur.def"
458                  stop
459               else
460                  nbq = nbq + 1
461                  niq(nbq) = i_cplus
462               end if
463            else
464               if (i_cplus /= 0) then
465                  write(*,*) "calchim: Error; C+ is present, but O2+ is not!!!"
466                  write(*,*) "Both must be in traceur.def if ionosphere wanted"
467                  stop
468               endif
469            endif   ! Of if (chemthermod == 3)
470            ! N+
471            i_nplus = igcm_nplus
472            if(chemthermod == 3) then
473               if (i_nplus == 0) then
474                  write(*,*) "calchim: Error; no N+ tracer !!!"
475                  write(*,*) "N+ is needed if O2+ is in traceur.def"
476                  stop
477               else
478                  nbq = nbq + 1
479                  niq(nbq) = i_nplus
480               end if
481            else
482               if (i_nplus /= 0) then
483                  write(*,*) "calchim: Error: N+ is present, but O2+ is not!!!"
484                  write(*,*) "Both must be in traceur.def if ionosphere wanted"
485                  stop
486               endif
487            endif   !Of if (chemthermod == 3)
488            ! NO+
489            i_noplus = igcm_noplus
490            if(chemthermod == 3) then
491               if (i_noplus == 0) then
492                  write(*,*) "calchim: Error; no NO+ tracer !!!"
493                  write(*,*) "NO+ is needed if O2+ is in traceur.def"
494                  stop
495               else
496                  nbq = nbq + 1
497                  niq(nbq) = i_noplus
498               end if
499            else
500               if (i_noplus /= 0) then
501                  write(*,*) "calchim: Error: NO+ is present, but O2+ is not!!!"
502                  write(*,*) "Both must be in traceur.def if ionosphere wanted"
503                  stop
504               endif
505            endif   !Of if (chemthermod == 3)
506            ! N2+
507            i_n2plus = igcm_n2plus
508            if (chemthermod == 3) then
509               if (i_n2plus == 0) then
510                  write(*,*) "calchim: Error; no N2+ tracer !!!"
511                  write(*,*) "N2+ is needed if O2+ is in traceur.def"
512                  stop
513               else
514                  nbq = nbq + 1
515                  niq(nbq) = i_n2plus
516               end if
517            else
518               if (i_n2plus /= 0) then
519                  write(*,*) "calchim: Error: N2+ is present, but O2+ is not!!!"
520                  write(*,*) "Both must be in traceur.def if ionosphere wanted"
521                  stop
522               endif
523            endif   !Of if (chemthermod == 3)
524            !H+
525            i_hplus = igcm_hplus
526            if (chemthermod == 3) then
527               if (i_hplus == 0) then
528                  write(*,*) "calchim: Error; no H+ tracer !!!"
529                  write(*,*) "H+ is needed if O2+ is in traceur.def"
530                  stop
531               else
532                  nbq = nbq + 1
533                  niq(nbq) = i_hplus
534               end if
535            else
536               if (i_hplus /= 0) then
537                  write(*,*) "calchim: Error: H+ is present, but O2+ is not!!!"
538                  write(*,*) "Both must be in traceur.def if ionosphere wanted"
539                  stop
540               endif
541            endif   !Of if (chemthermod == 3)
542            ! HCO2+
543            i_hco2plus = igcm_hco2plus
544            if(chemthermod == 3) then
545               if (i_hco2plus == 0) then
546                  write(*,*) "calchim: Error; no HCO2+ tracer !!!"
547                  write(*,*) "HCO2+ is needed if O2+ is in traceur.def"
548                  stop
549               else
550                  nbq = nbq + 1
551                  niq(nbq) = i_hco2plus
552               end if
553            else
554               if (i_hco2plus /= 0) then
555                  write(*,*) "calchim: Error: HCO2+ is present, but O2+ is not!!!"
556                  write(*,*) "Both must be in traceur.def if ionosphere wanted"
557                  stop
558               endif
559            endif    !Of if(chemthermod == 3)
560            !e-
561            i_elec = igcm_elec
562            if(chemthermod == 3) then
563               if (i_elec == 0) then
564                  write(*,*) "calchim: Error; no e- tracer !!!"
565                  write(*,*) "e- is needed if O2+ is in traceur.def"
566                  stop
567               else
568                  nbq = nbq + 1
569                  niq(nbq) = i_elec
570               end if
571            else
572               if(i_elec /= 0) then
573                  write(*,*) "calchim: Error: e- is present, but O2+ is not!!!"
574                  write(*,*) "Both must be in traceur.def if ionosphere wanted"
575                  stop
576               endif
577            endif   !Of if (chemthermod == 3)
578         endif      !Of thermochem
579
580         write(*,*) 'calchim: found nbq    = ',nbq,' tracers'
581               
582         firstcall = .false.
583      end if ! if (firstcall)
584
585! Initializations
586
587      zycol(:,:)    = 0.
588      dqchim(:,:,:) = 0.
589      dqschim(:,:)  = 0.
590
591!     latvl1= 22.27
592!     lonvl1= -47.94
593!     ig_vl1= 1+ int( (1.5-(latvl1-90.)*jjm/180.)  -2 )*iim +    &
594!             int(1.5+(lonvl1+180)*iim/360.)
595
596!=======================================================================
597!     loop over grid
598!=======================================================================
599
600      do ig = 1,ngrid
601         
602         foundswitch = 0
603         do l = 1,nlayer
604            do i = 1,nbq
605               iq = niq(i) ! get tracer index
606               zq(ig,l,iq) = pq(ig,l,iq) + pdq(ig,l,iq)*ptimestep
607               zycol(l,iq) = zq(ig,l,iq)*mmean(ig,l)/mmol(iq)
608            end do
609            zt(ig,l)  = pt(ig,l) + pdt(ig,l)*ptimestep
610            zu(ig,l)  = pu(ig,l) + pdu(ig,l)*ptimestep
611            zv(ig,l)  = pv(ig,l) + pdv(ig,l)*ptimestep
612            zpress(l) = pplay(ig,l)/100.
613            ztemp(l)  = zt(ig,l)
614            zdens(l)  = zpress(l)/(kb*1.e4*ztemp(l))
615            zlocal(l) = zzlay(ig,l)/1000.
616
617!           surfdust1d and surfice1d: conversion from m2/m3 to cm2/cm3
618
619            surfdust1d(l) = surfdust(ig,l)*1.e-2
620            surfice1d(l)  = surfice(ig,l)*1.e-2
621
622!           search for switch index between regions
623
624            if (photochem .and. thermochem) then
625               if (foundswitch == 0 .and. pplay(ig,l) < 1.e-1) then
626                  lswitch = l
627                  foundswitch = 1
628               end if
629            end if
630            if (.not. photochem) then
631               lswitch = 22
632            end if
633            if (.not. thermochem) then
634               lswitch = min(50,nlayer+1)
635            end if
636
637         end do ! of do l=1,nlayer
638
639         szacol = acos(mu0(ig))*180./pi
640         taucol = tauref(ig)*(700./610.)  ! provisoire en attente de nouveau jmars
641
642!=======================================================================
643!     call chemical subroutines
644!=======================================================================
645
646!        chemistry in lower atmosphere
647
648         if (photochem) then
649            call photochemistry(nlayer,nq,                         &
650                                lswitch,zycol,szacol,ptimestep,    &
651                                zpress,ztemp,zdens,dist_sol,       &
652                                surfdust1d,surfice1d,jo3,taucol)
653
654!        ozone photolysis, for output
655
656            do l = 1,nlayer
657               jo3_3d(ig,l) = jo3(l)
658            end do
659
660!        condensation of h2o2
661
662            call perosat(ngrid, nlayer, nq,                        &
663                         ig,ptimestep,pplev,pplay,                 &
664                         ztemp,zycol,dqcloud,dqscloud)
665         end if
666
667!        chemistry in upper atmosphere
668       
669         if (thermochem) then
670            call chemthermos(ig,lswitch,chemthermod,zycol,ztemp,zdens,  &
671                             zpress,zlocal,szacol,ptimestep,zday)
672         end if
673
674!        dry deposition
675
676         if (depos) then
677            call deposition(ngrid, nlayer, nq,                      &
678                            ig, ig_vl1, pplay, pplev, zzlay, zzlev, &
679                            zu, zv, zt, zycol, ptimestep, co2ice)
680         end if
681!=======================================================================
682!     tendencies
683!=======================================================================
684
685!     index of the most abundant species at each level
686
687!         major(:) = maxloc(zycol, dim = 2)
688
689!     tendency for the most abundant species = - sum of others
690         do l = 1,nlayer
691            iloc=maxloc(zycol(l,:))
692            iqmax=iloc(1)
693            do i = 1,nbq
694               iq = niq(i) ! get tracer index
695               if (iq /= iqmax) then
696                  dqchim(ig,l,iq) = (zycol(l,iq)*mmol(iq)/mmean(ig,l)  &
697                                   - zq(ig,l,iq))/ptimestep
698                  dqchim(ig,l,iqmax) = dqchim(ig,l,iqmax)              &
699                                     - dqchim(ig,l,iq)
700               end if
701            end do
702         end do ! of do l = 1,nlayer
703
704!=======================================================================
705!     end of loop over grid
706!=======================================================================
707
708      end do ! of do ig=1,ngrid
709
710!=======================================================================
711!     write outputs
712!=======================================================================
713
714! value of parameter 'output' to trigger writting of outputs
715! is set above at the declaration of the variable.
716
717      if (photochem .and. output) then
718         if (ngrid > 1) then
719            call writediagfi(ngrid,'jo3','j o3->o1d',    &
720                             's-1',3,jo3_3d(1,1))
721           if (callstats) then
722              call wstats(ngrid,'jo3','j o3->o1d',       &
723                          's-1',3,jo3_3d(1,1))
724           endif
725         end if ! of if (ngrid.gt.1)
726      end if ! of if (output)
727
728      return
729      end
730
Note: See TracBrowser for help on using the repository browser.