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

Last change on this file since 421 was 334, checked in by aslmd, 14 years ago

LMDZ.MARS:

28/10/11 == FL + AS

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