source: trunk/LMDZ.VENUS/libf/phyvenus/phytrac_chimie.F @ 3884

Last change on this file since 3884 was 3884, checked in by ikovalenko, 3 months ago
  • Property svn:executable set to *
File size: 18.7 KB
Line 
1      subroutine phytrac_chimie (
2     $                    debutphy,
3     $                    gmtime,
4     $                    nqmax,
5     $                    nlon,
6     $                    lat,
7     $                    lon,
8     $                    zzlay,
9     $                    nlev,
10     $                    pdtphys,
11     $                    temp,
12     $                    pplay,
13     $                    trac,
14     $                    d_tr_chem,
15     $                    iter,
16     $                    prod_tr,
17     $                    loss_tr,
18     $                    no_emi,
19     $                    o2_emi)
20
21      use chemparam_mod
22      use conc, only: mmean,rnew
23      use photolysis_mod, only: init_photolysis, nphot
24      use iono_h, only: temp_elect
25      use clesphys_mod
26      use YOMCST_mod
27#ifdef CPP_XIOS     
28      use xios_output_mod, only: send_xios_field
29#endif
30
31      implicit none
32     
33!#include "clesphys.h"
34!#include "YOMCST.h"
35
36!===================================================================
37!     input
38!===================================================================
39
40      integer :: nlon, nlev                     ! number of gridpoints and levels
41      integer :: nqmax                          ! number of tracers
42
43      real :: gmtime                            ! day fraction
44      real :: pdtphys                           ! phytrac_chimie timestep (s)
45      real, dimension(nlon,nlev) :: temp        ! temperature (k)
46      real, dimension(nlon,nlev) :: pplay       ! pressure (pa)
47      real, dimension(nlon,nlev) :: zzlay       ! altitude (m)
48      real, dimension(nlon,nlev,nqmax) :: trac  ! tracer mass mixing ratio
49
50      logical :: debutphy                       ! first call flag
51
52!===================================================================
53!     output
54!===================================================================
55
56      integer, dimension(nlon,nlev) :: iter                 ! chemical iterations
57      real, dimension(nlon,nlev,nqmax) :: d_tr_chem         ! chemical tendency for each tracer
58      real, dimension(nlon,nlev,nqmax) :: prod_tr, loss_tr  ! production and loss terms (for info)
59      real, dimension(nlon,nlev)    :: no_emi               ! no emission
60      real, dimension(nlon,nlev)    :: o2_emi               ! o2 emission
61
62!===================================================================
63!     local
64!===================================================================
65
66      real :: sza_local     ! solar zenith angle (deg)
67      real :: lon_sun
68      real :: zlocal(nlev)  ! altitude for photochem (km)
69      real :: t_elec(nlev)  ! electron temperature [K]
70     
71      integer, parameter :: t_elec_origin=2
72      !Electronic temperature. Index 1 -> Theis et al. 1980 - model data ; Index 2-> Theis et al. 1984 - model data
73     
74      integer :: i, iq
75      integer :: ilon, ilev
76
77      real  lat(nlon), lat_local(nlon)
78      real  lon(nlon), lon_local(nlon)
79
80      real, dimension(nlon,nlev) :: mrtwv, mrtsa ! total water and total sulfuric acid
81      real, dimension(nlon,nlev) :: mrwv, mrsa   ! gas-phase water and gas-phase sulfuric acid
82      real, dimension(nlon,nlev) :: trac_sum
83      real, dimension(nlon,nlev,nqmax) :: ztrac  ! local tracer mixing ratio
84      real, dimension(nlev) :: no_em
85      real, dimension(nlev) :: o2_em
86
87      integer, save :: nb_reaction_3_max     ! number of quadratic reactions
88      integer, save :: nb_reaction_4_max     ! number of bimolecular reactions
89      integer, save :: nquench               ! number of quenching + heterogeneous reactions
90      integer, save :: nphotion              ! number of photoionizations
91      integer, save :: nb_reaction_4_ion     ! quadratic reactions for ionosphere
92!      integer, save :: nb_reaction_4_deut    ! quadratic reactions for deuterium chem
93      integer, save :: nb_phot_max           ! total number of photolysis+photoionizations+quenching reactions
94
95      ! tracer indexes for the EUV heating:
96!!! ATTENTION. These values have to be identical to those in euvheat.F90
97!!! If the values are changed there, the same has to be done here  !!!
98
99      integer,parameter :: ix_co2  =  1
100      integer,parameter :: ix_co   =  2
101      integer,parameter :: ix_o    =  3
102      integer,parameter :: ix_o1d  =  4
103      integer,parameter :: ix_o2   =  5
104      integer,parameter :: ix_o3   =  6
105      integer,parameter :: ix_h    =  7
106      integer,parameter :: ix_h2   =  8
107      integer,parameter :: ix_oh   =  9
108      integer,parameter :: ix_ho2  = 10
109      integer,parameter :: ix_h2o2 = 11
110      integer,parameter :: ix_h2o  = 12
111      integer,parameter :: ix_n    = 13
112      integer,parameter :: ix_n2d  = 14
113      integer,parameter :: ix_no   = 15
114      integer,parameter :: ix_no2  = 16
115      integer,parameter :: ix_n2   = 17
116
117      ! NEED TO BE THE SAME AS IN EUVHEAT.F90
118      integer,parameter :: nespeuv = 17    ! Number of species considered (11, 12 or 17 (with nitrogen))
119
120      real :: vmr_dens_euv(nlon,nlev,nespeuv) ! local species density for EUV heating
121         
122!===================================================================
123!     first call : initialisations
124!===================================================================
125
126      if (debutphy) then
127     
128!--- Adjustment of Helium amount
129!       if (i_he/=0) then
130!          trac(:,:,i_he)=trac(:,:,i_he)/2.
131!       endif
132!---
133
134!-------------------------------------------------------------------
135!        Determination of chemistry reaction number
136!-------------------------------------------------------------------
137
138         if (ok_chem) then
139
140           ! set number of reactions, depending on ion chemistry or not
141
142           nb_reaction_4_ion  = 64
143           !nb_reaction_4_deut = 35
144   
145           ! default numbers if no ion and no deuterium chemistry
146   
147           nb_reaction_4_max = 99     ! set number of bimolecular reactions
148           nb_reaction_3_max = 14     ! set number of quadratic reactions
149           nquench           = 15     ! set number of first-order reactions:
150                                      !               quenching
151                                      !               thermal dissociation
152                                      !               heterogeneous
153           nphotion          = 0      ! set number of photoionizations
154   
155           if (ok_ionchem) then
156              nb_reaction_4_max = nb_reaction_4_max + nb_reaction_4_ion
157              nphotion       = 18     ! set number of photoionizations
158           endif
159
160           !if(deutchem) then
161           !   nb_reaction_4_max = nb_reaction_4_max + nb_reaction_4_deut 
162           !end if
163   
164           ! nb_phot_max is the total number of processes that are treated
165           ! numerically as a photolysis (first-order reaction):
166   
167           nb_phot_max = nphot + nphotion + nquench
168   
169!-------------------------------------------------------------------
170!        case of tracers re-initialisation with chemistry
171!-------------------------------------------------------------------
172 
173           if (reinit_trac .and. ok_chem) then
174 
175!             convert mass to volume mixing ratio
176              do iq = 1,nqmax - nmicro
177                 trac(:,:,iq) = trac(:,:,iq)*mmean(:,:)/m_tr(iq)
178              end do
179
180              print*, "SO2 is re-initialised"
181
182              if (i_so2 /= 0) then
183                 trac(:,25:,i_so2)  = 100.e-9
184                 trac(:,1:24,i_so2) = 10.e-6
185                 trac(:,25:,i_h2o)  = 1.e-6
186                 trac(:,1:24,i_h2o) = 30.e-6
187
188                 trac(:,:,i_osso_cis)   = 0.
189                 trac(:,:,i_osso_trans) = 0.
190                 trac(:,:,i_s2o2_cyc)   = 0.
191 
192! AL TRACERS!!!!!!!!!!!!!!!!!!!!!!!!!!!!
193!           print*, "Tracers are re-initialised"
194!           trac(:,:,:) = 1.0e-30
195 
196!             if ((i_ocs /= 0) .and. (i_co /= 0) .and. (i_hcl /= 0)
197!     $          .and. (i_so2 /= 0) .and. (i_h2o /= 0)
198!     $          .and. (i_n2 /= 0)  .and. (i_co2 /= 0)) then
199 
200!                trac(:,1:22,i_ocs) = 3.e-6
201!                trac(:,1:22,i_co)  = 25.e-6
202!                trac(:,:,i_hcl)    = 0.4e-6
203!                trac(:,1:22,i_so2) = 7.e-6
204!                trac(:,1:22,i_h2o) = 30.e-6
205!                trac(:,:,i_n2)     = 0.35e-1
206! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
207     
208!                ensure that sum of mixing ratios = 1
209 
210                 trac_sum(:,:) = 0.
211 
212                 do iq = 1,nqmax - nmicro
213                    if (iq /= i_co2) then
214                       trac_sum(:,:) = trac_sum(:,:) + trac(:,:,iq)
215                    end if
216                 end do
217 
218!                initialise co2
219 
220                 trac(:,:,i_co2) = 1. - trac_sum(:,:)
221 
222              else
223                 write(*,*) "at least one tracer is missing : stop"
224                 stop
225              end if
226         
227!             update mmean
228 
229              mmean(:,:) = 0.
230              do iq = 1,nqmax - nmicro
231                 mmean(:,:) = mmean(:,:)+trac(:,:,iq)*m_tr(iq)
232              enddo
233              rnew(:,:) = 8.314/mmean(:,:)*1.e3     ! J/kg K
234 
235!             convert volume to mass mixing ratio
236 
237              do iq = 1,nqmax - nmicro
238                 trac(:,:,iq) = trac(:,:,iq)*m_tr(iq)/mmean(:,:)
239              end do
240   
241           end if  ! reinit_trac
242
243         end if  ! ok_chem
244
245!-------------------------------------------------------------------
246!        case of detailed microphysics without chemistry
247!-------------------------------------------------------------------
248
249         if (.not. ok_chem .and. ok_cloud .and. cl_scheme == 2) then
250
251!           convert mass to volume mixing ratio
252
253            do iq = 1,nqmax - nmicro
254               ztrac(:,:,iq) = trac(:,:,iq)*mmean(:,:)/m_tr(iq)
255            end do
256         
257!           initialise microphysics
258 
259            call vapors4muphy_ini(nlon,nlev,ztrac)
260
261!           convert volume to mass mixing ratio
262
263            do iq = 1,nqmax - nmicro
264               trac(:,:,iq) = ztrac(:,:,iq)*m_tr(iq)/mmean(:,:)
265            end do
266   
267         end if
268
269      end if  ! debutphy
270
271!===================================================================
272!     convert mass to volume mixing ratio : gas phase
273!===================================================================
274
275      do iq = 1,nqmax - nmicro
276         ztrac(:,:,iq) = max(trac(:,:,iq)*mmean(:,:)/m_tr(iq), 1.e-30)
277      end do
278
279!===================================================================
280!     microphysics: simplified scheme (phd aurelien stolzenbach)
281!===================================================================
282
283      if (ok_cloud .and. cl_scheme == 1) then
284
285!     convert mass to volume mixing ratio : liquid phase
286
287         ztrac(:,:,i_h2so4liq) = max(trac(:,:,i_h2so4liq)
288     $                             *mmean(:,:)/m_tr(i_h2so4liq), 1.e-30)
289         ztrac(:,:,i_h2oliq) = max(trac(:,:,i_h2oliq)
290     $                             *mmean(:,:)/m_tr(i_h2oliq), 1.e-30)
291             
292!     total water and sulfuric acid (gas + liquid)
293
294         mrtwv(:,:) = ztrac(:,:,i_h2o) + ztrac(:,:,i_h2oliq)
295         mrtsa(:,:) = ztrac(:,:,i_h2so4) + ztrac(:,:,i_h2so4liq)
296
297!     all water and sulfuric acid is put in the gas-phase
298
299         mrwv(:,:) = mrtwv(:,:)
300         mrsa(:,:) = mrtsa(:,:)
301
302!     call microphysics
303
304         call new_cloud_venus(nb_mode, nlev, nlon, temp, pplay,
305     $                        mrtwv, mrtsa, mrwv, mrsa)
306
307!     update water vapour and sulfuric acid
308
309         ztrac(:,:,i_h2o) = mrwv(:,:)
310         ztrac(:,:,i_h2oliq) = mrtwv(:,:) - ztrac(:,:,i_h2o)
311     
312         ztrac(:,:,i_h2so4) = mrsa(:,:)
313         ztrac(:,:,i_h2so4liq) = mrtsa(:,:) - ztrac(:,:,i_h2so4)
314
315      end if  ! simplified scheme
316
317!===================================================================
318!     microphysics: detailed scheme (phd sabrina guilbon)
319!     !!! to be confirmed whether mad_muphy expects mmr or vmr for h2o and h2so4
320!===================================================================
321
322      if (ok_cloud .and. cl_scheme == 2) then
323
324         do iq = nqmax-nmicro+1,nqmax
325            ztrac(:,:,iq) = trac(:,:,iq)
326         end do
327
328         do ilon = 1,nlon       
329            do ilev = 1, nlev
330               if (temp(ilon,ilev) < 500.) then
331                  call mad_muphy(pdtphys,                               ! timestep
332     $                           temp(ilon,ilev),pplay(ilon,ilev),      ! temperature and pressure
333     $                           ztrac(ilon,ilev,i_h2o),
334     $                           ztrac(ilon,ilev,i_h2so4),     
335     $                           ztrac(ilon,ilev,i_m0_aer),
336     $                           ztrac(ilon,ilev,i_m3_aer),   
337     $                           ztrac(ilon,ilev,i_m0_mode1drop),
338     $                           ztrac(ilon,ilev,i_m0_mode1ccn),
339     $                           ztrac(ilon,ilev,i_m3_mode1sa),
340     $                           ztrac(ilon,ilev,i_m3_mode1w),   
341     $                           ztrac(ilon,ilev,i_m3_mode1ccn),   
342     $                           ztrac(ilon,ilev,i_m0_mode2drop),
343     $                           ztrac(ilon,ilev,i_m0_mode2ccn),
344     $                           ztrac(ilon,ilev,i_m3_mode2sa),
345     $                           ztrac(ilon,ilev,i_m3_mode2w),
346     $                           ztrac(ilon,ilev,i_m3_mode2ccn))
347               else
348                  ztrac(ilon,ilev,i_m0_aer)       = 0.
349                  ztrac(ilon,ilev,i_m3_aer)       = 0.
350                  ztrac(ilon,ilev,i_m0_mode1drop) = 0.
351                  ztrac(ilon,ilev,i_m0_mode1ccn)  = 0.
352                  ztrac(ilon,ilev,i_m3_mode1sa)   = 0.
353                  ztrac(ilon,ilev,i_m3_mode1w)    = 0.
354                  ztrac(ilon,ilev,i_m3_mode1ccn)  = 0.
355                  ztrac(ilon,ilev,i_m0_mode2drop) = 0.
356                  ztrac(ilon,ilev,i_m0_mode2ccn)  = 0.
357                  ztrac(ilon,ilev,i_m3_mode2sa)   = 0.
358                  ztrac(ilon,ilev,i_m3_mode2w)    = 0.
359                  ztrac(ilon,ilev,i_m3_mode2ccn)  = 0.
360               end if
361            end do
362         end do
363
364      end if  ! detailed scheme
365           
366!===================================================================
367!     photochemistry
368!===================================================================
369
370      if (ok_chem) then
371
372         lon_sun = (0.5 - gmtime)*2.*rpi
373         lon_local(:) = lon(:)*rpi/180.
374         lat_local(:) = lat(:)*rpi/180.
375         
376         if (ok_ionchem) then
377       
378           vmr_dens_euv(:,:,ix_co2) = ztrac(:,:,i_co2)  ! CO2
379           vmr_dens_euv(:,:,ix_co)  = ztrac(:,:,i_co)   ! CO
380           vmr_dens_euv(:,:,ix_o)   = ztrac(:,:,i_o)    ! O
381           vmr_dens_euv(:,:,ix_o1d) = ztrac(:,:,i_o1d)  ! O(1D)
382           vmr_dens_euv(:,:,ix_o2)  = ztrac(:,:,i_o2)   ! O2
383           vmr_dens_euv(:,:,ix_o3)  = ztrac(:,:,i_o3)   ! O3
384           vmr_dens_euv(:,:,ix_h)   = ztrac(:,:,i_h)    ! H
385           vmr_dens_euv(:,:,ix_h2)  = ztrac(:,:,i_h2)   ! H2
386           vmr_dens_euv(:,:,ix_oh)  = ztrac(:,:,i_oh)   ! OH
387           vmr_dens_euv(:,:,ix_ho2) = ztrac(:,:,i_ho2)  ! HO2
388           vmr_dens_euv(:,:,ix_h2o2)= ztrac(:,:,i_h2o2) ! H2O2
389           vmr_dens_euv(:,:,ix_h2o) = ztrac(:,:,i_h2o)  ! H2O
390           vmr_dens_euv(:,:,ix_n)   = ztrac(:,:,i_n)    ! N
391           vmr_dens_euv(:,:,ix_n2d) = ztrac(:,:,i_n2d)  ! N(2D)
392           vmr_dens_euv(:,:,ix_no)  = ztrac(:,:,i_no)   ! NO
393           vmr_dens_euv(:,:,ix_no2) = ztrac(:,:,i_no2)  ! NO2
394           vmr_dens_euv(:,:,ix_n2)  = ztrac(:,:,i_n2)   ! N2
395           
396         end if
397         
398         do ilon = 1,nlon
399            zlocal(:)=zzlay(ilon,:)/1000.
400                       
401!           solar zenith angle
402!            sza_local = acos(cos(lat_local(ilon))*cos(lon_local(ilon))
403!     $                 *cos(lon_sun) + cos(lat_local(ilon))
404!     $                 *sin(lon_local(ilon))*sin(lon_sun))*180./rpi
405
406            sza_local = cos(lat_local(ilon))*cos(lon_local(ilon))
407     $                 *cos(lon_sun) + cos(lat_local(ilon))
408     $                 *sin(lon_local(ilon))*sin(lon_sun)
409
410            ! Security - Handle rare cases where |sza_local| > 1
411            sza_local = min(sza_local,1.)
412            sza_local = max(-1.,sza_local)
413            sza_local = acos(sza_local)*180./rpi
414
415!           electron temperature
416
417            do ilev = 1, nlev
418               t_elec(ilev) = temp_elect(zlocal(ilev),temp(ilon,ilev),
419     $                                   sza_local,t_elec_origin) 
420            end do
421
422            call photochemistry_venus(nlev, nlon, zlocal, pdtphys,
423     $                           ok_jonline,ok_ionchem,tuneupperatm,
424     $                           nb_reaction_3_max,nb_reaction_4_max,
425     $                           nb_phot_max,nphotion,ilon,         
426     $                           pplay(ilon,:)/100.,
427     $                           temp(ilon,:), t_elec(:),
428     $                           ztrac(ilon,:,:),vmr_dens_euv(ilon,:,:),
429     $                           mmean(ilon,:),
430     $                           sza_local,
431     $                           lon(ilon), lat(ilon),
432     $                           nqmax, nespeuv, iter(ilon,:),
433     $                           prod_tr(ilon,:,:),
434     $                           loss_tr(ilon,:,:),
435     $                           no_em, o2_em)
436
437             no_emi(ilon,:) = no_em(:)
438             o2_emi(ilon,:) = o2_em(:)
439
440         end do         
441      end if  ! ok_chem
442
443!===================================================================
444!     compute tendencies
445!===================================================================
446
447!     update mmean
448
449      mmean(:,:) = 0.
450      do iq = 1,nqmax - nmicro
451         mmean(:,:) = mmean(:,:)+ztrac(:,:,iq)*m_tr(iq)
452      end do
453      rnew(:,:) = 8.314/mmean(:,:)*1.e3     ! J/kg K
454
455!===================================================================
456!     convert volume to mass mixing ratio / then tendencies in mmr
457!===================================================================
458
459!     gas phase
460
461      do iq = 1,nqmax - nmicro
462         ztrac(:,:,iq) = max(ztrac(:,:,iq)*m_tr(iq)/mmean(:,:),
463     $                       1.e-30)
464         d_tr_chem(:,:,iq) = (ztrac(:,:,iq) - trac(:,:,iq))/pdtphys
465      end do
466
467!     liquid phase or moments
468
469      if (ok_cloud .and. cl_scheme == 1) then
470         ztrac(:,:,i_h2so4liq) = max(ztrac(:,:,i_h2so4liq)
471     $                               *m_tr(i_h2so4liq)/mmean(:,:),
472     $                               1.e-30)
473         ztrac(:,:,i_h2oliq)   = max(ztrac(:,:,i_h2oliq)
474     $                               *m_tr(i_h2oliq)/mmean(:,:),
475     $                               1.e-30)
476         d_tr_chem(:,:,i_h2so4liq) = (ztrac(:,:,i_h2so4liq)
477     $                              - trac(:,:,i_h2so4liq))/pdtphys
478         d_tr_chem(:,:,i_h2oliq) = (ztrac(:,:,i_h2oliq)
479     $                            - trac(:,:,i_h2oliq))/pdtphys
480      end if
481
482
483      if (ok_cloud .and. cl_scheme == 2) then
484         do iq = nqmax-nmicro+1,nqmax
485            d_tr_chem(:,:,iq) = (ztrac(:,:,iq) - trac(:,:,iq))/pdtphys
486         end do
487      end if
488
489      end subroutine phytrac_chimie
Note: See TracBrowser for help on using the repository browser.