source: trunk/LMDZ.GENERIC/libf/aeronostd/calchim_asis.F90 @ 3567

Last change on this file since 3567 was 3312, checked in by yjaziri, 8 months ago

Generic PCM:

Photochemistry: correction for using diurnal equal true in 1D
No correction factor with diurnal equal true for photolysis rate
+ writediagspecUV change name to flux_surf for clarity

YJ

  • Property svn:executable set to *
File size: 21.1 KB
Line 
1      subroutine calchim_asis(ngrid,nlayer,nq,                              &
2                         ptimestep,pplay,pplev,pt,pdt,dist_sol,mu0,fract,   &
3                         zzlev,zzlay,zday,pq,pdq,dqchim,dqschim,            &
4                         tauref,co2ice,                                     &
5                         pu,pdu,pv,pdv,surfdust,surfice,icount,zdtchim)
6
7      use tracer_h,         only: mmol, noms, nesp, is_chim
8
9      use conc_mod,         only: mmean ! mean molecular mass of the atmosphere
10      USE comcstfi_mod
11      use callkeys_mod
12      use time_phylmdz_mod, only: ecritphy, iphysiq ! output rate set by ecritphy
13      use types_asis,       only: v_phot_3d, jlabel, nb_phot_hv_max
14      use chimiedata_h
15      use radcommon_h,      only: glat
16      use wstats_mod,       only: wstats
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!    update 20/10/2017 adapted to LMDZ GENERIC+cosmetic changes (Benjamin Charnay)
38!    update 06/03/2021 generic tracer/network + photolysis online (Yassin Jaziri)
39!
40!   Arguments:
41!   ----------
42!
43!  Input:
44!
45!    ptimestep                  timestep (s)
46!    pplay(ngrid,nlayer)    Pressure at the middle of the layers (Pa)
47!    pplev(ngrid,nlayer+1)  Intermediate pressure levels (Pa)
48!    pt(ngrid,nlayer)       Temperature (K)
49!    pdt(ngrid,nlayer)      Temperature tendency (K)
50!    pu(ngrid,nlayer)       u component of the wind (ms-1)
51!    pdu(ngrid,nlayer)      u component tendency (K)
52!    pv(ngrid,nlayer)       v component of the wind (ms-1)
53!    pdv(ngrid,nlayer)      v component tendency (K)
54!    dist_sol                   distance of the sun (AU)
55!    mu0(ngrid)               cos of solar zenith angle (=1 when sun at zenith)
56!    fract(ngrid)           day fraction (for diurnal=false)
57!    pq(ngrid,nlayer,nq)  Advected fields, ie chemical species here
58!    pdq(ngrid,nlayer,nq) Previous tendencies on pq
59!    tauref(ngrid)            Optical depth at 7 hPa
60!    co2ice(ngrid)            co2 ice surface layer (kg.m-2)
61!    surfdust(ngrid,nlayer) dust surface area (m2/m3)
62!    surfice(ngrid,nlayer)  ice surface area (m2/m3)
63!
64!  Output:
65!
66!    dqchim(ngrid,nlayer,nesp)   ! tendencies on pq due to chemistry
67!    dqschim(ngrid,nesp)         ! tendencies on qsurf
68!
69!=======================================================================
70
71!     input:
72
73      integer,intent(in) :: ngrid    ! number of atmospheric columns
74      integer,intent(in) :: nlayer   ! number of atmospheric layers
75      integer,intent(in) :: nq       ! number of tracers
76      real :: ptimestep
77      real :: pplay(ngrid,nlayer)    ! pressure at the middle of the layers
78      real :: zzlay(ngrid,nlayer)    ! altitude at the middle of the layers
79      real :: pplev(ngrid,nlayer+1)  ! intermediate pressure levels
80      real :: zzlev(ngrid,nlayer+1)  ! altitude at layer boundaries
81      real :: pt(ngrid,nlayer)       ! temperature
82      real :: pdt(ngrid,nlayer)      ! temperature tendency
83      real :: pu(ngrid,nlayer)       ! u component of the wind (m.s-1)
84      real :: pdu(ngrid,nlayer)      ! u component tendency
85      real :: pv(ngrid,nlayer)       ! v component of the wind (m.s-1)
86      real :: pdv(ngrid,nlayer)      ! v component tendency
87      real :: dist_sol               ! distance of the sun (AU)
88      real :: mu0(ngrid)             ! cos of solar zenith angle (=1 when sun at zenith)
89      real :: pq(ngrid,nlayer,nq)    ! tracers mass mixing ratio
90      real :: pdq(ngrid,nlayer,nq)   ! previous tendencies
91      real :: zday                   ! date (time since Ls=0, in martian days)
92      real :: tauref(ngrid)          ! optical depth at 7 hPa
93      real :: co2ice(ngrid)          ! co2 ice surface layer (kg.m-2)
94      real :: surfdust(ngrid,nlayer) ! dust surface area (m2/m3)
95      real :: surfice(ngrid,nlayer)  !  ice surface area (m2/m3)
96
97!     output:
98
99      real :: dqchim(ngrid,nlayer,nesp)   ! tendencies on pq due to chemistry
100      real :: dqschim(ngrid,nesp)         ! tendencies on qsurf
101
102!     local variables:
103
104      integer :: iloc(1)                  ! index of major species
105      integer :: ig,l,i,iq,iqmax,iesp
106      integer :: foundswitch, lswitch     ! to switch between photochem and thermochem ? (YJ)
107      integer,save :: chemthermod
108
109      real    :: zq(ngrid,nlayer,nesp) ! pq+pdq*ptimestep before chemistry
110                                       ! new mole fraction after
111      real    :: zt(ngrid,nlayer)      ! temperature
112      real    :: zu(ngrid,nlayer)      ! u component of the wind
113      real    :: zv(ngrid,nlayer)      ! v component of the wind
114      real    :: taucol                ! optical depth at 7 hPa
115      real    :: xmmol(nlayer,nesp)    ! mmol/mmean but only for chemical species
116 
117      logical,save :: firstcall = .true.
118
119!     for each column of atmosphere:
120
121      real :: zpress(nlayer)       ! Pressure (mbar)
122      real :: zdens(nlayer)        ! Density  (cm-3)
123      real :: ztemp(nlayer)        ! Temperature (K)
124      real :: zlocal(nlayer)       ! Altitude (km)
125      real :: zycol(nlayer,nesp)   ! Composition (mole fractions)
126      real :: zmmean(nlayer)       ! Mean molecular mass (g.mole-1)
127      real :: szacol               ! Solar zenith angle
128      real :: fract(ngrid)         ! day fraction
129      real :: fractcol             ! day fraction
130      real :: surfice1d(nlayer)    ! Ice surface area (cm2/cm3)
131      real :: surfdust1d(nlayer)   ! Dust surface area (cm2/cm3)
132
133      integer :: iter(nlayer)      !  Number of chemical iterations
134                                   !  within one physical timestep
135      integer :: icount
136!     for output:
137
138      logical :: output              ! to issue calls to writediagfi and stats
139      parameter (output = .true.)
140      real :: iter_3d(ngrid,nlayer)  ! Number of chemical iterations
141                                     !  within one physical timestep
142      integer           ::  ierr
143      real              ::  zdtchim(ngrid,nlayer)    ! temperature modification by chemistry
144      real              ::  dEzchim(ngrid,nlayer)    ! energy modification by chemistry
145      real              ::  zdtchim_output(ngrid)    ! flux modification by chemistry in W.m-2
146
147!=======================================================================
148!     initialization of the chemistry (first call only)
149!=======================================================================
150
151      if (firstcall) then
152
153!! Moved to routine indice in photochemistry_asis
154!! because nb_phot_hv_max value needed in order
155!! to choose if we call read_phototable or not.
156!! A cleaner solution need to be find.
157!         if (photochem .and. .not. jonline) then
158!           print*,'calchim: Read photolysis lookup table'
159!           call read_phototable
160!         end if
161
162         if (.not.allocated(SF_mode))       allocate(SF_mode(nesp))
163         if (.not.allocated(SF_value))      allocate(SF_value(nesp))
164         if (.not.allocated(prod_rate))     allocate(prod_rate(nesp))
165         if (.not.allocated(surface_flux))  allocate(surface_flux(ngrid,nesp))
166         if (.not.allocated(surface_flux2)) allocate(surface_flux2(ngrid,nesp))
167         if (.not.allocated(escape))        allocate(escape(ngrid,nesp))
168         if (.not.allocated(chemnoms))      allocate(chemnoms(nesp))
169
170         surface_flux(:,:)  = 0.
171         surface_flux2(:,:) = 0.
172         escape(:,:)        = 0.
173         SF_mode(:)         = 2
174         SF_value(:)        = 0.
175         prod_rate(:)       = 0.
176         iter_3d(:,:)       = 0.
177         iter(:)            = 0.
178
179         call ini_tracchim
180
181         ! Sanity check mmol /= 0. in chemistry
182         do iq = 1,nq
183            if (is_chim(iq).eq.1 .and. mmol(iq).eq.0.) then
184               write(*,*) 'Error in calchim:'
185               write(*,*) 'Mmol cannot be equal to 0 for chemical species'
186               stop
187            end if
188         end do
189
190         firstcall = .false.
191      end if ! if (firstcall)
192
193! Initializations
194
195      zycol(:,:)    = 0.
196      dqchim(:,:,:) = 0.
197      dqschim(:,:)  = 0.
198
199!=======================================================================
200!     loop over grid
201!=======================================================================
202
203      do ig = 1,ngrid
204         
205         foundswitch = 0
206         do l = 1,nlayer
207            iesp = 0
208            do iq = 1,nq
209               if (is_chim(iq).eq.1) then
210                  iesp = iesp + 1
211                  zq(ig,l,iesp) = pq(ig,l,iq) + pdq(ig,l,iq)*ptimestep
212                  xmmol(l,iesp) = mmol(iq)/mmean(ig,l)
213                  zycol(l,iesp) = zq(ig,l,iesp)/xmmol(l,iesp)
214               end if
215            end do
216
217            zt(ig,l)  = pt(ig,l) + pdt(ig,l)*ptimestep
218            zu(ig,l)  = pu(ig,l) + pdu(ig,l)*ptimestep
219            zv(ig,l)  = pv(ig,l) + pdv(ig,l)*ptimestep
220            zpress(l) = pplay(ig,l)/100.
221            ztemp(l)  = zt(ig,l)
222            zdens(l)  = zpress(l)/(kb*1.e4*ztemp(l))
223            zlocal(l) = zzlay(ig,l)/1000.
224            zmmean(l) = mmean(ig,l)
225
226!           surfdust1d and surfice1d: conversion from m2/m3 to cm2/cm3
227
228            surfdust1d(l) = surfdust(ig,l)*1.e-2
229            surfice1d(l)  = surfice(ig,l)*1.e-2
230
231!           search for switch index between regions
232
233!            if (photochem .and. thermochem) then
234!               if (foundswitch == 0 .and. pplay(ig,l) < 1.e-1) then
235!                  lswitch = l
236!                  foundswitch = 1
237!               end if
238!            end if
239            if (.not. photochem) then
240               lswitch = 22
241            end if
242!            if (.not. thermochem) then
243               lswitch = min(50,nlayer+1)
244!            end if
245
246         end do ! of do l=1,nlayer
247
248         szacol   = acos(mu0(ig))*180./pi
249         taucol   = tauref(ig)*(700./610.)  ! provisoire en attente de nouveau jmars
250         fractcol = fract(ig)
251
252!=======================================================================
253!     call chemical subroutines
254!=======================================================================
255
256!        chemistry in lower atmosphere
257
258            call photochemistry_asis(nlayer,ngrid,                         &
259                                ig,lswitch,zycol,szacol,fractcol,ptimestep,&
260                                zpress,zlocal,ztemp,zdens,zmmean,dist_sol, &
261                                surfdust1d,surfice1d,taucol,iter,zdtchim(ig,:))
262
263            ! diagnostic photochemical heating
264            zdtchim_output(ig) = 0.
265            do l = 1,nlayer
266              dEzchim(ig,l) = zdtchim(ig,l)*cpp*(pplev(ig,l) - pplev(ig,l+1))/glat(ig)
267              zdtchim_output(ig) = zdtchim_output(ig) + zdtchim(ig,l)*cpp*(pplev(ig,l) - pplev(ig,l+1))/glat(ig)
268            end do
269
270            do l = 1,nlayer
271               iter_3d(ig,l) = iter(l)
272            end do
273
274!        chemistry in upper atmosphere
275       
276!         if (thermochem) then
277!            call chemthermos(ig,nlayer,lswitch,chemthermod,zycol,ztemp,&
278!                             zdens,zpress,zlocal,szacol,ptimestep,zday)
279!         end if
280
281!        dry deposition
282
283         if (depos) then
284            call deposition_source(ngrid, nlayer, nq,   &
285                        ig, zzlay, zzlev, zdens, zycol, ptimestep)
286            surface_flux2(ig,:) = surface_flux2(ig,:) + surface_flux(ig,:)
287            if (ngrid==1) then
288              if(mod(icount,ecritphy).eq.0) then
289                surface_flux2(ig,:) = surface_flux2(ig,:)/float(ecritphy)
290              endif
291            else
292              if(mod(icount*iphysiq,ecritphy).eq.0) then
293                surface_flux2(ig,:) = surface_flux2(ig,:)*iphysiq/float(ecritphy)
294              endif
295            endif
296         end if
297
298!=======================================================================
299!     tendencies
300!=======================================================================
301
302!     index of the most abundant species at each level
303
304!         major(:) = maxloc(zycol, dim = 2)
305
306!     tendency for the most abundant species = - sum of others
307
308         do l = 1,nlayer
309            iloc=maxloc(zycol(l,:))
310            iqmax=iloc(1)
311            do iq = 1,nesp
312               if (iq /= iqmax) then
313                  dqchim(ig,l,iq) = (zycol(l,iq)*xmmol(l,iq) - zq(ig,l,iq))/ptimestep
314                  dqchim(ig,l,iq) = max(dqchim(ig,l,iq),-zq(ig,l,iq)/ptimestep)
315                  dqchim(ig,l,iqmax) = dqchim(ig,l,iqmax) - dqchim(ig,l,iq)
316               end if
317            end do
318
319         end do ! of do l = 1,nlayer
320
321
322
323
324
325
326!=======================================================================
327!     end of loop over grid
328!=======================================================================
329
330      end do ! of do ig=1,ngridbidon(ig,:) = 1
331
332!=======================================================================
333!     write outputs
334!=======================================================================
335
336! value of parameter 'output' to trigger writting of outputs
337! is set above at the declaration of the variable.
338      if (output) then
339
340            ! photoloysis
341            do i=1,nb_phot_hv_max
342               call wstats(ngrid,trim(jlabel(i,1)),jlabel(i,1),    &
343                                's-1',3,v_phot_3d(1,1,i))
344            end do
345            ! iter
346            call wstats(ngrid,'iter','iterations',  &
347                             ' ',3,iter_3d(1,1))
348            ! phothcemical heating
349            call wstats(ngrid,'zdtchim','dT chim',    &
350                             'K.s-1',3,zdtchim)
351            call wstats(ngrid,'dEzchim','dE chim',    &
352                             'W m-2',3,dEzchim)
353            call wstats(ngrid,"Ezchim","integrated dT chim","W m-2",2,zdtchim_output)
354            ! Mean molecular mass
355            call wstats(ngrid,'mmean','mean molecular mass',       &
356                           'g.mole-1',3,mmean(1,1))
357            ! Chemical tendencies
358            do iesp=1,nesp
359               call wstats(ngrid,trim(chemnoms(iesp))//'_dq',trim(chemnoms(iesp))//'_dq',    &
360                               'kg/kg s^-1',3,dqchim(1,1,iesp) )
361            end do
362            if (depos) then
363               ! Surface fluxes and escape
364               do iesp=1,nesp
365                  call wstats(ngrid,trim(chemnoms(iesp))//'_flux_mean',trim(chemnoms(iesp))//' mean flux',    &
366                               'molecule.m-2.s-1',2,surface_flux2(1,indexchim(trim(chemnoms(iesp)))))
367                  call wstats(ngrid,trim(chemnoms(iesp))//'_flux',trim(chemnoms(iesp))//' flux',    &
368                               'molecule.m-2.s-1',2,surface_flux(1,indexchim(trim(chemnoms(iesp)))))
369                  call wstats(ngrid,trim(chemnoms(iesp))//'_escape',trim(chemnoms(iesp))//' escape',    &
370                               'molecule.m-2.s-1',2,escape(1,indexchim(trim(chemnoms(iesp)))))
371               end do
372            endif ! end depos
373
374         ! photoloysis
375         do i=1,nb_phot_hv_max
376            call writediagfi(ngrid,trim(jlabel(i,1)),jlabel(i,1),    &
377                             's-1',3,v_phot_3d(1,1,i))
378         end do
379         ! iter
380         call writediagfi(ngrid,'iter','iterations',  &
381                          ' ',3,iter_3d(1,1))
382         ! phothcemical heating
383         call writediagfi(ngrid,'zdtchim','dT chim',    &
384                          'K.s-1',3,zdtchim)
385         call writediagfi(ngrid,'dEzchim','dE chim',    &
386                          'W m-2',3,dEzchim)
387         call writediagfi(ngrid,"Ezchim","integrated dT chim","W m-2",2,zdtchim_output)
388         ! Mean molecular mass
389         call writediagfi(ngrid,'mmean','mean molecular mass',       &
390                        'g.mole-1',3,mmean(1,1))
391         ! Chemical tendencies
392         do iesp=1,nesp
393            call writediagfi(ngrid,trim(chemnoms(iesp))//'_dq',trim(chemnoms(iesp))//'_dq',    &
394                            'kg/kg s^-1',3,dqchim(1,1,iesp) )
395         end do
396         if (depos) then
397            ! Surface fluxes and escape
398            do iesp=1,nesp
399               call writediagfi(ngrid,trim(chemnoms(iesp))//'_flux_mean',trim(chemnoms(iesp))//' mean flux',    &
400                            'molecule.m-2.s-1',2,surface_flux2(1,indexchim(trim(chemnoms(iesp)))))
401               call writediagfi(ngrid,trim(chemnoms(iesp))//'_flux',trim(chemnoms(iesp))//' flux',    &
402                            'molecule.m-2.s-1',2,surface_flux(1,indexchim(trim(chemnoms(iesp)))))
403               call writediagfi(ngrid,trim(chemnoms(iesp))//'_escape',trim(chemnoms(iesp))//' escape',    &
404                            'molecule.m-2.s-1',2,escape(1,indexchim(trim(chemnoms(iesp)))))
405            end do
406            ! Restart flux average
407            if (ngrid == 1) then
408               if(mod(icount,ecritphy).eq.0) then
409                  surface_flux2(:,:) = 0.0
410               endif
411            else
412               if(mod(icount*iphysiq,ecritphy).eq.0) then
413                  surface_flux2(:,:) = 0.0
414               endif
415            endif
416         endif ! end depos
417
418         if (jonline .and. output_writediagspecUV) then
419
420            !flux at the surface
421            call writediagspecUV(ngrid,"flux_surf","Surface flux(lon,lat,band)","photon.s-1.nm-1.cm-2",3,fluxUV(:,:,1))
422
423         endif ! end jonline
424
425      end if ! of if (output)
426
427      return
428
429
430      contains
431
432
433!======================================================================
434
435      subroutine ini_tracchim
436
437!======================================================================
438! YJ Modern tracer.def
439! Get in tracer.def chemical parameters
440!======================================================================
441
442         use chimiedata_h
443         use tracer_h, only:   is_chim
444         implicit none
445
446         ! local
447         character(len=500) :: tracline ! to store a line of text
448         integer            :: nq       ! number of tracers
449         integer            :: iesp, iq
450
451         ! Get nq
452         open(90,file='traceur.def',status='old',form='formatted',iostat=ierr)
453         if (ierr.eq.0) then
454           READ(90,'(A)') tracline
455           IF (trim(tracline).ne.'#ModernTrac-v1') THEN ! Test modern traceur.def
456            READ(tracline,*,iostat=ierr) nq ! Try standard traceur.def
457           ELSE
458              DO
459                READ(90,'(A)',iostat=ierr) tracline
460                IF (ierr.eq.0) THEN
461                  IF (index(tracline,'#').ne.1) THEN ! Allows arbitary number of comments lines in the header
462                    READ(tracline,*,iostat=ierr) nq
463                    EXIT
464                  ENDIF
465                ELSE ! If pb, or if reached EOF without having found number of tracer
466                   write(*,*) "calchim: error reading line of tracers"
467                   write(*,*) "   (first lines of traceur.def) "
468                   stop
469                ENDIF
470              ENDDO
471           ENDIF ! if modern or standard traceur.def
472         else
473            write(*,*) "calchim: error opening traceur.def"
474            stop
475         endif
476
477         ! Get data of tracers
478         iesp = 0
479         do iq=1,nq
480            read(90,'(A)') tracline
481            if (is_chim(iq).eq.1) then
482               iesp = iesp + 1
483               read(tracline,*) chemnoms(iesp)
484               write(*,*)"ini_tracchim: iq=",iq,"noms(iq)=",trim(chemnoms(iesp))
485               if (index(tracline,'SF_mode='    ) /= 0) then
486                  read(tracline(index(tracline,'SF_mode=')+len('SF_mode='):),*) SF_mode(iesp)
487                  write(*,*) ' Parameter value (traceur.def) : SF_mode=', SF_mode(iesp)
488               else
489                  write(*,*) ' Parameter value (default)     : SF_mode=', SF_mode(iesp)
490               end if
491               if (index(tracline,'SF_value='   ) /= 0) then
492                  read(tracline(index(tracline,'SF_value=')+len('SF_value='):),*) SF_value(iesp)
493                  write(*,*) ' Parameter value (traceur.def) : SF_value=', SF_value(iesp)
494               else
495                  write(*,*) ' Parameter value (default)     : SF_value=', SF_value(iesp)
496               end if
497               if (index(tracline,'prod_rate='  ) /= 0) then
498                  read(tracline(index(tracline,'prod_rate=')+len('prod_rate='):),*) prod_rate(iesp)
499                  write(*,*) ' Parameter value (traceur.def) : prod_rate=', prod_rate(iesp)
500               else
501                  write(*,*) ' Parameter value (default)     : prod_rate=', prod_rate(iesp)
502               end if
503               if (index(tracline,'escape='     ) /= 0) then
504                  read(tracline(index(tracline,'escape=')+len('escape='):),*) escape(:,iesp)
505                  write(*,*) ' Parameter value (traceur.def) : escape=', escape(:,iesp)
506               else
507                  write(*,*) ' Parameter value (default)     : escape=', escape(:,iesp)
508               end if
509            end if
510         end do
511
512         close(90)
513
514      end subroutine ini_tracchim
515
516      end subroutine calchim_asis
517
Note: See TracBrowser for help on using the repository browser.