source: trunk/LMDZ.MARS/libf/aeronomars/calchim.F @ 553

Last change on this file since 553 was 476, checked in by emillour, 14 years ago

Mars GCM: minor correction in aeronomars/calchim.F and aeronomars/surfacearea.F : calls to wstats should always be embeded in if (callstats) condition.
EM

File size: 16.8 KB
Line 
1      subroutine calchim(ptimestep,pplay,pplev,pt,pdt,dist_sol,mu0,
2     $                   zzlev,zzlay,zday,pq,pdq,dqchim,dqschim,dqcloud,
3     $                   dqscloud,tauref,co2ice,
4     $                   pu,pdu,pv,pdv,surfdust,surfice)
5c
6      implicit none
7c
8c=======================================================================
9c
10c   subject:
11c   --------
12c
13c  Prepare the call for the photochemical module, and send back the
14c  tendencies from photochemistry in the chemical species mass mixing ratios
15c
16c   Author:   Sebastien Lebonnois (08/11/2002)
17c   -------
18c    update 12/06/2003 for water ice clouds and compatibility with dust
19c    update 07/2003 for coupling with thermosphere (Monica Angelats-i-Coll)
20c    update 03/05/2005 cosmetic changes (Franck Lefevre)
21c    update sept. 2008 identify tracers by their names (Ehouarn Millour)
22c    update 17/03/2011 synchronize with latest version of chemistry (Franck Lefevre)
23c    update 05/12/2011 synchronize with latest version of chemistry (Franck Lefevre)
24c
25c   Arguments:
26c   ----------
27c
28c  Input:
29c
30c    ptimestep                  timestep (s)
31c    pplay(ngridmx,nlayermx)    Pressure at the middle of the layers (Pa)
32c    pplev(ngridmx,nlayermx+1)  Intermediate pressure levels (Pa)
33c    pt(ngridmx,nlayermx)       Temperature (K)
34c    pdt(ngridmx,nlayermx)      Temperature tendency (K)
35c    pu(ngridmx,nlayermx)       u component of the wind (ms-1)
36c    pdu(ngridmx,nlayermx)      u component tendency (K)
37c    pv(ngridmx,nlayermx)       v component of the wind (ms-1)
38c    pdv(ngridmx,nlayermx)      v component tendency (K)
39c    dist_sol                   distance of the sun (AU)
40c    mu0(ngridmx)               cos of solar zenith angle (=1 when sun at zenith)
41c    pq(ngridmx,nlayermx,nqmx)  Advected fields, ie chemical species here
42c    pdq(ngridmx,nlayermx,nqmx) Previous tendencies on pq
43c    tauref(ngridmx)            Optical depth at 7 hPa
44c    co2ice(ngridmx)            co2 ice surface layer (kg.m-2)
45c    surfdust(ngridmx,nlayermx) dust surface area (m2/m3)
46c    surfice(ngridmx,nlayermx)  ice surface area (m2/m3)
47c
48c  Output:
49c
50c    dqchim(ngridmx,nlayermx,nqmx) ! tendencies on pq due to chemistry
51c    dqschim(ngridmx,nqmx)         ! tendencies on qsurf
52c
53c=======================================================================
54
55c    Declarations :
56c    --------------
57
58#include "dimensions.h"
59#include "dimphys.h"
60#include "chimiedata.h"
61#include "tracer.h"
62#include "comcstfi.h"
63#include "callkeys.h"
64#include "conc.h"
65
66c Arguments :
67c -----------
68
69c   inputs:
70c   -------
71
72      real    ptimestep
73      real    pplay(ngridmx,nlayermx)    ! pressure at the middle of the layers
74      real    zzlay(ngridmx,nlayermx)    ! pressure at the middle of the layers
75      real    pplev(ngridmx,nlayermx+1)  ! intermediate pressure levels
76      real    zzlev(ngridmx,nlayermx+1)  ! altitude at layer boundaries
77      real    pt(ngridmx,nlayermx)       ! temperature
78      real    pdt(ngridmx,nlayermx)      ! temperature tendency
79      real    pu(ngridmx,nlayermx)       ! u component of the wind (m.s-1)
80      real    pdu(ngridmx,nlayermx)      ! u component tendency
81      real    pv(ngridmx,nlayermx)       ! v component of the wind (m.s-1)
82      real    pdv(ngridmx,nlayermx)      ! v component tendency
83      real    dist_sol                   ! distance of the sun (AU)
84      real    mu0(ngridmx)               ! cos of solar zenith angle (=1 when sun at zenith)
85      real    pq(ngridmx,nlayermx,nqmx)  ! tracers mass mixing ratio
86      real    pdq(ngridmx,nlayermx,nqmx) ! previous tendencies
87      real    zday                       ! date (time since Ls=0, in martian days)
88      real    tauref(ngridmx)            ! optical depth at 7 hPa
89      real    co2ice(ngridmx)            ! co2 ice surface layer (kg.m-2)
90      real    surfdust(ngridmx,nlayermx) ! dust surface area (m2/m3)
91      real    surfice(ngridmx,nlayermx)  !  ice surface area (m2/m3)
92
93c   outputs:
94c   --------
95
96      real dqchim(ngridmx,nlayermx,nqmx) ! tendencies on pq due to chemistry
97      real dqschim(ngridmx,nqmx)         ! tendencies on qsurf
98      real dqcloud(ngridmx,nlayermx,nqmx)! tendencies on pq due to condensation
99      real dqscloud(ngridmx,nqmx)        ! tendencies on qsurf
100
101c Local variables :
102c -----------------
103
104      character*20 str20
105      integer  ig,l,i,iq
106      integer  foundswitch, lswitch
107      integer  ig_vl1
108      real latvl1, lonvl1
109      real zq(ngridmx,nlayermx,nqmx) ! pq+pdq*ptimestep before chemistry
110                                     ! new mole fraction after
111      real zt(ngridmx,nlayermx)      ! temperature
112      real zu(ngridmx,nlayermx)      ! u component of the wind
113      real zv(ngridmx,nlayermx)      ! v component of the wind
114      real taucol                    ! optical depth at 7 hPa
115
116      logical depos                  ! switch for dry deposition
117      parameter (depos=.false.)
118c
119c  for each column of atmosphere:
120c
121      real zpress(nlayermx)       !  Pressure (mbar)
122      real zdens(nlayermx)        !  Density  (cm-3)
123      real ztemp(nlayermx)        !  Temperature (K)
124      real zlocal(nlayermx)       !  Altitude (km)
125      real zycol(nlayermx,nqmx)   !  Composition (mole fractions)
126      real szacol                 !  Solar zenith angle
127      real surfice1d(nlayermx)    !  Ice surface area (cm2/cm3)
128      real surfdust1d(nlayermx)   !  Dust surface area (cm2/cm3)
129      real jo3(nlayermx)          !  Photodissociation rate O3->O1D (s-1)
130c
131c for output:
132c
133      real jo3_3d(ngridmx,nlayermx)  !  Photodissociation rate O3->O1D (s-1)
134
135      logical output              ! to issue calls to writediagfi and stats
136      parameter (output=.true.)   ! see at end of routine
137
138      logical,save :: firstcall=.true.
139      integer,save :: nbq       ! number of tracers used in the chemistry
140      integer,save :: niq(nqmx) ! array storing the indexes of the tracers
141
142! index of tracers:
143
144      integer,save :: i_co2=0
145      integer,save :: i_co=0
146      integer,save :: i_o=0
147      integer,save :: i_o1d=0
148      integer,save :: i_o2=0
149      integer,save :: i_o3=0
150      integer,save :: i_h=0
151      integer,save :: i_h2=0
152      integer,save :: i_oh=0
153      integer,save :: i_ho2=0
154      integer,save :: i_h2o2=0
155      integer,save :: i_ch4=0
156      integer,save :: i_n2=0
157      integer,save :: i_ar=0
158      integer,save :: i_ice=0 ! water ice
159      integer,save :: i_h2o=0 ! water vapour
160c
161c=======================================================================
162c     initialization of the chemistry (first call only)
163c=======================================================================
164c
165      if (firstcall) then
166c
167         if (photochem) then
168            print*,'calchim: Read photolysis lookup table'
169            call read_phototable
170         end if
171
172         ! find index of chemical tracers to use
173         nbq=0 ! to count number of tracers
174         i_co2=igcm_co2
175         if (i_co2.eq.0) then
176           write(*,*) "calchim: Error; no CO2 tracer !!!"
177           stop
178         else
179           nbq=nbq+1
180           niq(nbq)=i_co2
181         endif
182         i_co=igcm_co
183         if (i_co.eq.0) then
184           write(*,*) "calchim: Error; no CO tracer !!!"
185           stop
186         else
187           nbq=nbq+1
188           niq(nbq)=i_co
189         endif
190         i_o=igcm_o
191         if (i_o.eq.0) then
192           write(*,*) "calchim: Error; no O tracer !!!"
193           stop
194         else
195           nbq=nbq+1
196           niq(nbq)=i_o
197         endif
198         i_o1d=igcm_o1d
199         if (i_o1d.eq.0) then
200           write(*,*) "calchim: Error; no O1D tracer !!!"
201           stop
202         else
203           nbq=nbq+1
204           niq(nbq)=i_o1d
205         endif
206         i_o2=igcm_o2
207         if (i_o2.eq.0) then
208           write(*,*) "calchim: Error; no O2 tracer !!!"
209           stop
210         else
211           nbq=nbq+1
212           niq(nbq)=i_o2
213         endif
214         i_o3=igcm_o3
215         if (i_o3.eq.0) then
216           write(*,*) "calchim: Error; no O3 tracer !!!"
217           stop
218         else
219           nbq=nbq+1
220           niq(nbq)=i_o3
221         endif
222         i_h=igcm_h
223         if (i_h.eq.0) then
224           write(*,*) "calchim: Error; no H tracer !!!"
225           stop
226         else
227           nbq=nbq+1
228           niq(nbq)=i_h
229         endif
230         i_h2=igcm_h2
231         if (i_h2.eq.0) then
232           write(*,*) "calchim: Error; no H2 tracer !!!"
233           stop
234         else
235           nbq=nbq+1
236           niq(nbq)=i_h2
237         endif
238         i_oh=igcm_oh
239         if (i_oh.eq.0) then
240           write(*,*) "calchim: Error; no OH tracer !!!"
241           stop
242         else
243           nbq=nbq+1
244           niq(nbq)=i_oh
245         endif
246         i_ho2=igcm_ho2
247         if (i_ho2.eq.0) then
248           write(*,*) "calchim: Error; no HO2 tracer !!!"
249           stop
250         else
251           nbq=nbq+1
252           niq(nbq)=i_ho2
253         endif
254         i_h2o2=igcm_h2o2
255         if (i_h2o2.eq.0) then
256           write(*,*) "calchim: Error; no H2O2 tracer !!!"
257           stop
258         else
259           nbq=nbq+1
260           niq(nbq)=i_h2o2
261         endif
262         i_ch4=igcm_ch4
263         if (i_ch4.eq.0) then
264           write(*,*) "calchim: Error; no CH4 tracer !!!"
265           stop
266         else
267           nbq=nbq+1
268           niq(nbq)=i_ch4
269         endif
270         i_n2=igcm_n2
271         if (i_n2.eq.0) then
272           write(*,*) "calchim: Error; no N2 tracer !!!"
273           stop
274         else
275           nbq=nbq+1
276           niq(nbq)=i_n2
277         endif
278         i_ar=igcm_ar
279         if (i_ar.eq.0) then
280           write(*,*) "calchim: Error; no AR tracer !!!"
281           stop
282         else
283           nbq=nbq+1
284           niq(nbq)=i_ar
285         endif
286         i_ice=igcm_h2o_ice
287         if (i_ice.eq.0) then
288           write(*,*) "calchim: Error; no water ice tracer !!!"
289           stop
290         else
291           nbq=nbq+1
292           niq(nbq)=i_ice
293         endif
294         i_h2o=igcm_h2o_vap
295         if (i_h2o.eq.0) then
296           write(*,*) "calchim: Error; no water vapor tracer !!!"
297           stop
298         else
299           nbq=nbq+1
300           niq(nbq)=i_h2o
301         endif
302
303         write(*,*) 'calchim: found nbq    = ',nbq,' tracers'
304         write(*,*) '               i_co2  = ',i_co2
305         write(*,*) '               i_co   = ',i_co
306         write(*,*) '               i_o    = ',i_o
307         write(*,*) '               i_o1d  = ',i_o1d
308         write(*,*) '               i_o2   = ',i_o2
309         write(*,*) '               i_o3   = ',i_o3
310         write(*,*) '               i_h    = ',i_h
311         write(*,*) '               i_h2   = ',i_h2
312         write(*,*) '               i_oh   = ',i_oh
313         write(*,*) '               i_ho2  = ',i_ho2
314         write(*,*) '               i_h2o2 = ',i_h2o2
315         write(*,*) '               i_ch4  = ',i_ch4
316         write(*,*) '               i_n2   = ',i_n2
317         write(*,*) '               i_ar   = ',i_ar
318         write(*,*) '               i_ice  = ',i_ice
319         write(*,*) '               i_h2o  = ',i_h2o
320         
321         firstcall = .false.
322      end if ! if (firstcall)
323
324! Initialize output tendencies to zero (to handle case of tracers which
325! are not used in the chemistry (e.g. dust))
326
327      dqchim(:,:,:) = 0
328      dqschim(:,:)  = 0
329c
330!     latvl1= 22.27
331!     lonvl1= -47.94
332!     ig_vl1= 1+ int( (1.5-(latvl1-90.)*jjm/180.)  -2 )*iim +
333!    $        int(1.5+(lonvl1+180)*iim/360.)
334c
335c=======================================================================
336c     loop over grid
337c=======================================================================
338c
339      do ig = 1,ngridmx
340c
341c     local updates
342c
343         foundswitch = 0
344         do l = 1,nlayermx
345            do i = 1,nbq
346              iq = niq(i) ! get tracer index
347              zq(ig,l,iq) = pq(ig,l,iq) + pdq(ig,l,iq)*ptimestep
348              zycol(l,iq) = zq(ig,l,iq)*mmean(ig,l)/mmol(iq)
349            end do
350            zt(ig,l) = pt(ig,l) + pdt(ig,l)*ptimestep
351            zu(ig,l) = pu(ig,l) + pdu(ig,l)*ptimestep
352            zv(ig,l) = pv(ig,l) + pdv(ig,l)*ptimestep
353c
354            zpress(l) = pplay(ig,l)/100.
355            ztemp(l)  = zt(ig,l)
356            zdens(l)  = zpress(l)/(kb*1.e4*ztemp(l))
357            zlocal(l) = zzlay(ig,l)/1000.
358c
359c           surfdust1d and surfice1d: conversion from m2/m3 to cm2/cm3
360c
361            surfdust1d(l) = surfdust(ig,l)*1.e-2
362            surfice1d(l)  = surfice(ig,l)*1.e-2
363c
364c           search for switch index between regions
365c
366            if (photochem .and. thermochem) then
367               if (foundswitch .eq. 0 .and. pplay(ig,l).lt.1.e-3) then
368                  lswitch = l
369                  foundswitch=1
370               end if
371            end if
372            if ( .not. photochem) then
373               lswitch = 22
374            end if
375            if (.not.  thermochem) then
376               lswitch = min(50,nlayermx+1)    ! FL (original value: 33)
377            end if
378c
379         end do ! of do l=1,nlayermx
380c
381         szacol = acos(mu0(ig))*180./pi
382         taucol = tauref(ig)
383c
384c=======================================================================
385c     call chemical subroutine
386c=======================================================================
387c
388c        chemistry in lower atmosphere
389c
390         if (photochem) then
391            call photochemistry(lswitch,zycol,szacol,ptimestep,
392     $                          zpress,ztemp,zdens,dist_sol,
393     $                          surfdust1d,surfice1d,jo3,taucol)
394         end if
395c
396c        chemistry in upper atmosphere
397c
398         if (thermochem) then
399            call chemthermos(ig,lswitch,zycol,ztemp,zdens,zpress,
400     $                       zlocal,szacol,ptimestep,zday)
401         end if
402c
403c        dry deposition
404c
405         if (depos) then
406            call deposition(ig, ig_vl1, pplay, pplev, zzlay, zzlev,
407     $                      zu, zv, zt, zycol, ptimestep, co2ice)
408         end if
409c
410c=======================================================================
411c     tendencies
412c=======================================================================
413c
414c     must be 0. for water ice:
415c
416         if (water) then
417            do l = 1,nlayermx
418               dqchim(ig,l,i_ice) = 0.
419            end do
420         end if
421c
422c     tendency for CO2 = - sum of others for lower atmosphere
423c     tendency for O   = - sum of others for upper atmosphere
424c
425         do l = 1,nlayermx
426            if (l .lt. lswitch) then
427               do i=1,nbq
428                 iq=niq(i) ! get tracer index
429                  if ((iq.ne.i_co2).and.(iq.ne.i_ice)) then
430                     dqchim(ig,l,iq) = (zycol(l,iq)*mmol(iq)/mmean(ig,l)
431     $                                - zq(ig,l,iq))/ptimestep
432                  else if (iq.eq.i_co2) then
433                     dqchim(ig,l,iq) = 0.
434                  end if
435                  dqschim(ig,iq) = 0.
436               end do ! of do i=1,nbq
437               
438               do i=1,nbq
439                 iq=niq(i) ! get tracer index
440                  if (iq.ne.i_co2) then
441                     dqchim(ig,l,i_co2) = dqchim(ig,l,i_co2)
442     $                                  - dqchim(ig,l,iq)
443                  end if
444               end do
445             else if (l .ge. lswitch) then
446               do i=1,nbq
447                 iq=niq(i) ! get tracer index
448                   if ((iq.ne.i_o).and.(iq.ne.i_ice)) then
449                      dqchim(ig,l,iq) = (zycol(l,iq)*mmol(iq)
450     $                                  /mmean(ig,l)
451     $                                 - zq(ig,l,iq))/ptimestep
452                   else if (iq.eq.i_o) then
453                      dqchim(ig,l,iq) = 0.
454                   end if
455               enddo
456               
457               do i=1,nbq
458                 iq=niq(i) ! get tracer index
459                   if (iq.ne.i_o) then
460                      dqchim(ig,l,i_o) = dqchim(ig,l,i_o)
461     $                                 - dqchim(ig,l,iq)
462                   end if
463               end do
464             end if ! of if (l.lt.lswitch) else if (l.ge.lswitch)
465          end do ! of do l = 1,nlayermx
466c
467c     condensation of h2o2
468c
469          call perosat(ig,ptimestep,pplev,pplay,
470     $                 ztemp,zycol,dqcloud,dqscloud)
471c
472c     density and j(o3->o1d), for outputs
473c
474          do l = 1,nlayermx
475             jo3_3d(ig,l) = jo3(l)
476          end do
477c
478c=======================================================================
479c     end of loop over grid
480c=======================================================================
481c
482      end do ! of do ig=1,ngridmx
483c
484c=======================================================================
485c     write outputs
486c=======================================================================
487c
488! value of parameter 'output' to trigger writting of outputs
489! is set above at the declaration of the variable.
490
491      if (output) then
492         if (ngridmx .gt. 1) then
493            call writediagfi(ngridmx,'jo3','j o3->o1d',
494     $                       's-1',3,jo3_3d(1,1))
495           if (callstats) then
496            call wstats(ngridmx,'jo3','j o3->o1d',
497     $                       's-1',3,jo3_3d(1,1))
498           endif
499         end if ! of if (ngridmx.gt.1)
500      endif ! of if (output)
501c
502      end
Note: See TracBrowser for help on using the repository browser.