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

Last change on this file since 3981 was 3893, checked in by gmilcareck, 4 months ago

Remove all "call abort" and "stop" statement in LMDZ.GENERIC and replacing them by call abort_physic().

  • Property svn:executable set to *
File size: 20.8 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: diagfi_output_rate ! output rate
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,intent(in) :: icount ! physics iteration number
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               call abort_physic("calchim","Mmol cannot be equal to 0 for chemical species",1)
185            end if
186         end do
187
188         firstcall = .false.
189      end if ! if (firstcall)
190
191! Initializations
192
193      zycol(:,:)    = 0.
194      dqchim(:,:,:) = 0.
195      dqschim(:,:)  = 0.
196
197!=======================================================================
198!     loop over grid
199!=======================================================================
200
201      do ig = 1,ngrid
202         
203         foundswitch = 0
204         do l = 1,nlayer
205            iesp = 0
206            do iq = 1,nq
207               if (is_chim(iq).eq.1) then
208                  iesp = iesp + 1
209                  zq(ig,l,iesp) = pq(ig,l,iq) + pdq(ig,l,iq)*ptimestep
210                  xmmol(l,iesp) = mmol(iq)/mmean(ig,l)
211                  zycol(l,iesp) = zq(ig,l,iesp)/xmmol(l,iesp)
212               end if
213            end do
214
215            zt(ig,l)  = pt(ig,l) + pdt(ig,l)*ptimestep
216            zu(ig,l)  = pu(ig,l) + pdu(ig,l)*ptimestep
217            zv(ig,l)  = pv(ig,l) + pdv(ig,l)*ptimestep
218            zpress(l) = pplay(ig,l)/100.
219            ztemp(l)  = zt(ig,l)
220            zdens(l)  = zpress(l)/(kb*1.e4*ztemp(l))
221            zlocal(l) = zzlay(ig,l)/1000.
222            zmmean(l) = mmean(ig,l)
223
224!           surfdust1d and surfice1d: conversion from m2/m3 to cm2/cm3
225
226            surfdust1d(l) = surfdust(ig,l)*1.e-2
227            surfice1d(l)  = surfice(ig,l)*1.e-2
228
229!           search for switch index between regions
230
231!            if (photochem .and. thermochem) then
232!               if (foundswitch == 0 .and. pplay(ig,l) < 1.e-1) then
233!                  lswitch = l
234!                  foundswitch = 1
235!               end if
236!            end if
237            if (.not. photochem) then
238               lswitch = 22
239            end if
240!            if (.not. thermochem) then
241               lswitch = min(50,nlayer+1)
242!            end if
243
244         end do ! of do l=1,nlayer
245
246         szacol   = acos(mu0(ig))*180./pi
247         taucol   = tauref(ig)*(700./610.)  ! provisoire en attente de nouveau jmars
248         fractcol = fract(ig)
249
250!=======================================================================
251!     call chemical subroutines
252!=======================================================================
253
254!        chemistry in lower atmosphere
255
256            call photochemistry_asis(nlayer,ngrid,                         &
257                                ig,lswitch,zycol,szacol,fractcol,ptimestep,&
258                                zpress,zlocal,ztemp,zdens,zmmean,dist_sol, &
259                                surfdust1d,surfice1d,taucol,iter,zdtchim(ig,:))
260
261            ! diagnostic photochemical heating
262            zdtchim_output(ig) = 0.
263            do l = 1,nlayer
264              dEzchim(ig,l) = zdtchim(ig,l)*cpp*(pplev(ig,l) - pplev(ig,l+1))/glat(ig)
265              zdtchim_output(ig) = zdtchim_output(ig) + zdtchim(ig,l)*cpp*(pplev(ig,l) - pplev(ig,l+1))/glat(ig)
266            end do
267
268            do l = 1,nlayer
269               iter_3d(ig,l) = iter(l)
270            end do
271
272!        chemistry in upper atmosphere
273       
274!         if (thermochem) then
275!            call chemthermos(ig,nlayer,lswitch,chemthermod,zycol,ztemp,&
276!                             zdens,zpress,zlocal,szacol,ptimestep,zday)
277!         end if
278
279!        dry deposition
280
281         if (depos) then
282            call deposition_source(ngrid, nlayer, nq,   &
283                        ig, zzlay, zzlev, zdens, zycol, ptimestep)
284            surface_flux2(ig,:) = surface_flux2(ig,:) + surface_flux(ig,:)
285            if(mod(icount,diagfi_output_rate).eq.0) then
286              !note that surface_flux2(:,:) is reset every diagfi_output_rate
287              surface_flux2(ig,:) = surface_flux2(ig,:)/float(diagfi_output_rate)
288            endif
289         end if ! of if (depos)
290
291!=======================================================================
292!     tendencies
293!=======================================================================
294
295!     index of the most abundant species at each level
296
297!         major(:) = maxloc(zycol, dim = 2)
298
299!     tendency for the most abundant species = - sum of others
300
301         do l = 1,nlayer
302            iloc=maxloc(zycol(l,:))
303            iqmax=iloc(1)
304            do iq = 1,nesp
305               if (iq /= iqmax) then
306                  dqchim(ig,l,iq) = (zycol(l,iq)*xmmol(l,iq) - zq(ig,l,iq))/ptimestep
307                  dqchim(ig,l,iq) = max(dqchim(ig,l,iq),-zq(ig,l,iq)/ptimestep)
308                  dqchim(ig,l,iqmax) = dqchim(ig,l,iqmax) - dqchim(ig,l,iq)
309               end if
310            end do
311
312         end do ! of do l = 1,nlayer
313
314
315
316
317
318
319!=======================================================================
320!     end of loop over grid
321!=======================================================================
322
323      end do ! of do ig=1,ngridbidon(ig,:) = 1
324
325!=======================================================================
326!     write outputs
327!=======================================================================
328
329! value of parameter 'output' to trigger writting of outputs
330! is set above at the declaration of the variable.
331      if (output) then
332
333            ! photoloysis
334            do i=1,nb_phot_hv_max
335               call wstats(ngrid,trim(jlabel(i,1)),jlabel(i,1),    &
336                                's-1',3,v_phot_3d(1,1,i))
337            end do
338            ! iter
339            call wstats(ngrid,'iter','iterations',  &
340                             ' ',3,iter_3d(1,1))
341            ! phothcemical heating
342            call wstats(ngrid,'zdtchim','dT chim',    &
343                             'K.s-1',3,zdtchim)
344            call wstats(ngrid,'dEzchim','dE chim',    &
345                             'W m-2',3,dEzchim)
346            call wstats(ngrid,"Ezchim","integrated dT chim","W m-2",2,zdtchim_output)
347            ! Mean molecular mass
348            call wstats(ngrid,'mmean','mean molecular mass',       &
349                           'g.mole-1',3,mmean(1,1))
350            ! Chemical tendencies
351            do iesp=1,nesp
352               call wstats(ngrid,trim(chemnoms(iesp))//'_dq',trim(chemnoms(iesp))//'_dq',    &
353                               'kg/kg s^-1',3,dqchim(1,1,iesp) )
354            end do
355            if (depos) then
356               ! Surface fluxes and escape
357               do iesp=1,nesp
358                  call wstats(ngrid,trim(chemnoms(iesp))//'_flux_mean',trim(chemnoms(iesp))//' mean flux',    &
359                               'molecule.m-2.s-1',2,surface_flux2(1,indexchim(trim(chemnoms(iesp)))))
360                  call wstats(ngrid,trim(chemnoms(iesp))//'_flux',trim(chemnoms(iesp))//' flux',    &
361                               'molecule.m-2.s-1',2,surface_flux(1,indexchim(trim(chemnoms(iesp)))))
362                  call wstats(ngrid,trim(chemnoms(iesp))//'_escape',trim(chemnoms(iesp))//' escape',    &
363                               'molecule.m-2.s-1',2,escape(1,indexchim(trim(chemnoms(iesp)))))
364               end do
365            endif ! end depos
366
367         ! photoloysis
368         do i=1,nb_phot_hv_max
369            call writediagfi(ngrid,trim(jlabel(i,1)),jlabel(i,1),    &
370                             's-1',3,v_phot_3d(1,1,i))
371         end do
372         ! iter
373         call writediagfi(ngrid,'iter','iterations',  &
374                          ' ',3,iter_3d(1,1))
375         ! phothcemical heating
376         call writediagfi(ngrid,'zdtchim','dT chim',    &
377                          'K.s-1',3,zdtchim)
378         call writediagfi(ngrid,'dEzchim','dE chim',    &
379                          'W m-2',3,dEzchim)
380         call writediagfi(ngrid,"Ezchim","integrated dT chim","W m-2",2,zdtchim_output)
381         ! Mean molecular mass
382         call writediagfi(ngrid,'mmean','mean molecular mass',       &
383                        'g.mole-1',3,mmean(1,1))
384         ! Chemical tendencies
385         do iesp=1,nesp
386            call writediagfi(ngrid,trim(chemnoms(iesp))//'_dq',trim(chemnoms(iesp))//'_dq',    &
387                            'kg/kg s^-1',3,dqchim(1,1,iesp) )
388         end do
389         if (depos) then
390            ! Surface fluxes and escape
391            do iesp=1,nesp
392               call writediagfi(ngrid,trim(chemnoms(iesp))//'_flux_mean',trim(chemnoms(iesp))//' mean flux',    &
393                            'molecule.m-2.s-1',2,surface_flux2(1,indexchim(trim(chemnoms(iesp)))))
394               call writediagfi(ngrid,trim(chemnoms(iesp))//'_flux',trim(chemnoms(iesp))//' flux',    &
395                            'molecule.m-2.s-1',2,surface_flux(1,indexchim(trim(chemnoms(iesp)))))
396               call writediagfi(ngrid,trim(chemnoms(iesp))//'_escape',trim(chemnoms(iesp))//' escape',    &
397                            'molecule.m-2.s-1',2,escape(1,indexchim(trim(chemnoms(iesp)))))
398            end do
399            ! Restart flux average
400            if(mod(icount,diagfi_output_rate).eq.0) then
401              surface_flux2(:,:) = 0.0
402            endif
403         endif ! end depos
404
405         if (jonline .and. output_writediagspecUV) then
406
407            !flux at the surface
408            call writediagspecUV(ngrid,"flux_surf","Surface flux(lon,lat,band)","photon.s-1.nm-1.cm-2",3,fluxUV(:,:,1))
409
410         endif ! end jonline
411
412      end if ! of if (output)
413
414
415      contains
416
417
418!======================================================================
419
420      subroutine ini_tracchim
421
422!======================================================================
423! YJ Modern tracer.def
424! Get in tracer.def chemical parameters
425!======================================================================
426
427         use chimiedata_h
428         use tracer_h, only:   is_chim
429         implicit none
430
431         ! local
432         character(len=500) :: tracline ! to store a line of text
433         integer            :: nq       ! number of tracers
434         integer            :: iesp, iq
435
436         ! Get nq
437         open(90,file='traceur.def',status='old',form='formatted',iostat=ierr)
438         if (ierr.eq.0) then
439           READ(90,'(A)') tracline
440           IF (trim(tracline).ne.'#ModernTrac-v1') THEN ! Test modern traceur.def
441            READ(tracline,*,iostat=ierr) nq ! Try standard traceur.def
442           ELSE
443              DO
444                READ(90,'(A)',iostat=ierr) tracline
445                IF (ierr.eq.0) THEN
446                  IF (index(tracline,'#').ne.1) THEN ! Allows arbitary number of comments lines in the header
447                    READ(tracline,*,iostat=ierr) nq
448                    EXIT
449                  ENDIF
450                ELSE ! If pb, or if reached EOF without having found number of tracer
451                   write(*,*) "calchim: error reading line of tracers"
452                   write(*,*) "   (first lines of traceur.def) "
453                   call abort_physic("calchim", "error reading line of tracers",1)
454                ENDIF
455              ENDDO
456           ENDIF ! if modern or standard traceur.def
457         else
458            call abort_physic("calchim", "Unable to open file traceur.def",1)
459         endif
460
461         ! Get data of tracers
462         iesp = 0
463         do iq=1,nq
464            read(90,'(A)') tracline
465            if (is_chim(iq).eq.1) then
466               iesp = iesp + 1
467               read(tracline,*) chemnoms(iesp)
468               write(*,*)"ini_tracchim: iq=",iq,"noms(iq)=",trim(chemnoms(iesp))
469               if (index(tracline,'SF_mode='    ) /= 0) then
470                  read(tracline(index(tracline,'SF_mode=')+len('SF_mode='):),*) SF_mode(iesp)
471                  write(*,*) ' Parameter value (traceur.def) : SF_mode=', SF_mode(iesp)
472               else
473                  write(*,*) ' Parameter value (default)     : SF_mode=', SF_mode(iesp)
474               end if
475               if (index(tracline,'SF_value='   ) /= 0) then
476                  read(tracline(index(tracline,'SF_value=')+len('SF_value='):),*) SF_value(iesp)
477                  write(*,*) ' Parameter value (traceur.def) : SF_value=', SF_value(iesp)
478               else
479                  write(*,*) ' Parameter value (default)     : SF_value=', SF_value(iesp)
480               end if
481               if (index(tracline,'prod_rate='  ) /= 0) then
482                  read(tracline(index(tracline,'prod_rate=')+len('prod_rate='):),*) prod_rate(iesp)
483                  write(*,*) ' Parameter value (traceur.def) : prod_rate=', prod_rate(iesp)
484               else
485                  write(*,*) ' Parameter value (default)     : prod_rate=', prod_rate(iesp)
486               end if
487               if (index(tracline,'escape='     ) /= 0) then
488                  read(tracline(index(tracline,'escape=')+len('escape='):),*) escape(:,iesp)
489                  write(*,*) ' Parameter value (traceur.def) : escape=', escape(:,iesp)
490               else
491                  write(*,*) ' Parameter value (default)     : escape=', escape(:,iesp)
492               end if
493            end if
494         end do
495
496         close(90)
497
498      end subroutine ini_tracchim
499
500      end subroutine calchim_asis
501
Note: See TracBrowser for help on using the repository browser.