source: trunk/LMDZ.MARS/libf/aeronomars/euvheat.F90 @ 2042

Last change on this file since 2042 was 2042, checked in by emillour, 6 years ago

Mars GCM:
Modifications to use the parametrized photoabsorbtion coefficients;
a first step towards implementing ionospheric chemistry in the new
chemical solver:

  • change in species indexes in chemthermos.F90, paramfoto_compact.F, hrtherm.F and euvheat.F90
  • calchim.F90: added a variable in call to photochemistry
  • photochemistry.F90: added calls to jthermcalc_e107 and phdisrate, with an additionlal flag, jparam (.false. by default). The computed photodissociation coefficents are sent to v_phot, which is used in the chemistry. Thus concentrations computed in chimtogcm are now done over all atmospheric layers.

FGG

File size: 14.7 KB
Line 
1      SUBROUTINE euvheat(ngrid,nlayer,nq,pt,pdt,pplev,pplay,zzlay, &
2           mu0,ptimestep,ptime,zday,pq,pdq,pdteuv)
3
4      use tracer_mod, only: igcm_co2, igcm_co, igcm_o, igcm_o1d,         &
5                            igcm_o2, igcm_h, igcm_h2, igcm_oh, igcm_ho2, &
6                            igcm_h2o2, igcm_h2o_vap, igcm_o3, igcm_n2,   &
7                            igcm_n, igcm_no, igcm_no2, igcm_n2d, mmol
8      use conc_mod, only: rnew, cpnew
9      IMPLICIT NONE
10!=======================================================================
11!   subject:
12!   --------
13!   Computing heating rate due to EUV absorption
14!
15!   author:  MAC 2002
16!   ------
17!
18!   input:
19!   -----
20!   mu0(ngrid)           
21!   pplay(ngrid,nlayer)   pressure at middle of layers (Pa)
22!
23!   output:
24!   -------
25!
26!   pdteuv(ngrid,nlayer)      Heating rate (K/s)
27!
28!=======================================================================
29!
30!    0.  Declarations :
31!    ------------------
32!
33#include "callkeys.h"
34!-----------------------------------------------------------------------
35!    Input/Output
36!    ------------
37
38      integer,intent(in) :: ngrid ! number of atmospheric columns
39      integer,intent(in) :: nlayer ! number of atmospheric layers
40      integer,intent(in) :: nq ! number of advected tracers
41      real :: pt(ngrid,nlayer)
42      real :: pdt(ngrid,nlayer)
43      real :: pplev(ngrid,nlayer+1)
44      real :: pplay(ngrid,nlayer)
45      real :: zzlay(ngrid,nlayer)
46      real :: mu0(ngrid)
47      real :: ptimestep,ptime
48      real :: zday
49      real :: pq(ngrid,nlayer,nq)
50      real :: pdq(ngrid,nlayer,nq)
51
52      real :: pdteuv(ngrid,nlayer)
53!
54!    Local variables :
55!    -----------------
56      integer,save :: nespeuv    ! Number of species considered
57
58      INTEGER :: l,ig,n
59      integer,save :: euvmod
60      real, allocatable, save :: rm(:,:)   ! number density (cm-3)
61      real :: zq(ngrid,nlayer,nq) ! local updated tracer quantity
62      real :: zt(ngrid,nlayer)      ! local updated atmospheric temperature
63      real :: zlocal(nlayer)
64      real :: zenit
65      real :: jtot(nlayer)
66      real :: dens                                              ! amu/cm-3
67      real :: tx(nlayer)
68!      real euveff     !UV heating efficiency
69     
70! tracer indexes for the EUV heating:
71!!! ATTENTION. These values have to be identical to those in chemthermos.F90
72!!! If the values are changed there, the same has to be done here  !!!
73      integer,parameter :: i_co2  =  1
74      integer,parameter :: i_co   =  2
75      integer,parameter :: i_o    =  3
76      integer,parameter :: i_o1d  =  4
77      integer,parameter :: i_o2   =  5
78      integer,parameter :: i_o3   =  6
79      integer,parameter :: i_h    =  7
80      integer,parameter :: i_h2   =  8
81      integer,parameter :: i_oh   =  9
82      integer,parameter :: i_ho2  = 10
83      integer,parameter :: i_h2o2 = 11
84      integer,parameter :: i_h2o  = 12
85      integer,parameter :: i_n    = 13
86      integer,parameter :: i_n2d  = 14
87      integer,parameter :: i_no   = 15
88      integer,parameter :: i_no2  = 16
89      integer,parameter :: i_n2   = 17
90!      integer,parameter :: i_co2=1
91!      integer,parameter :: i_o2=2
92!      integer,parameter :: i_o=3
93!      integer,parameter :: i_co=4
94!      integer,parameter :: i_h=5
95!      integer,parameter :: i_oh=6
96!      integer,parameter :: i_ho2=7
97!      integer,parameter :: i_h2=8
98!      integer,parameter :: i_h2o=9
99!      integer,parameter :: i_h2o2=10
100!      integer,parameter :: i_o1d=11
101!      integer,parameter :: i_o3=12
102!      integer,parameter :: i_n2=13
103!      integer,parameter :: i_n=14
104!      integer,parameter :: i_no=15
105!      integer,parameter :: i_n2d=16
106!      integer,parameter :: i_no2=17
107
108     
109! Tracer indexes in the GCM:
110      integer,save :: g_co2=0
111      integer,save :: g_o=0
112      integer,save :: g_o2=0
113      integer,save :: g_h2=0
114      integer,save :: g_h2o2=0
115      integer,save :: g_h2o=0
116      integer,save :: g_o3=0
117      integer,save :: g_n2=0
118      integer,save :: g_n=0
119      integer,save :: g_no=0
120      integer,save :: g_co=0
121      integer,save :: g_h=0
122      integer,save :: g_no2=0
123      integer,save :: g_oh=0
124      integer,save :: g_ho2=0
125      integer,save :: g_o1d=0
126      integer,save :: g_n2d=0
127
128
129      logical,save :: firstcall=.true.
130
131! Initializations and sanity checks:
132
133
134      if (firstcall) then
135         nespeuv=0
136        ! identify the indexes of the tracers we'll need
137         g_co2=igcm_co2
138         if (g_co2.eq.0) then
139            write(*,*) "euvheat: Error; no CO2 tracer !!!"
140            write(*,*) "CO2 is always needed if calleuv=.true."
141            stop
142         else
143            nespeuv=nespeuv+1
144         endif
145         g_o=igcm_o
146         if (g_o.eq.0) then
147            write(*,*) "euvheat: Error; no O tracer !!!"
148            write(*,*) "O is always needed if calleuv=.true."
149            stop
150         else
151            nespeuv=nespeuv+1
152         endif
153         g_o2=igcm_o2
154         if (g_o2.eq.0) then
155            write(*,*) "euvheat: Error; no O2 tracer !!!"
156            write(*,*) "O2 is always needed if calleuv=.true."
157            stop
158         else
159            nespeuv=nespeuv+1
160         endif
161         g_h2=igcm_h2
162         if (g_h2.eq.0) then
163            write(*,*) "euvheat: Error; no H2 tracer !!!"
164            write(*,*) "H2 is always needed if calleuv=.true."
165            stop
166         else
167            nespeuv=nespeuv+1
168         endif
169         g_oh=igcm_oh
170         if (g_oh.eq.0) then
171            write(*,*) "euvheat: Error; no OH tracer !!!"
172            write(*,*) "OH must always be present if thermochem=T"
173            stop
174         else
175            nespeuv=nespeuv+1 
176         endif
177         g_ho2=igcm_ho2
178         if (g_ho2.eq.0) then
179            write(*,*) "euvheat: Error; no HO2 tracer !!!"
180            write(*,*) "HO2 must always be present if thermochem=T"
181            stop
182         else
183            nespeuv=nespeuv+1 
184         endif
185         g_h2o2=igcm_h2o2
186         if (g_h2o2.eq.0) then
187            write(*,*) "euvheat: Error; no H2O2 tracer !!!"
188            write(*,*) "H2O2 is always needed if calleuv=.true."
189            stop
190         else
191            nespeuv=nespeuv+1
192         endif
193         g_h2o=igcm_h2o_vap
194         if (g_h2o.eq.0) then
195            write(*,*) "euvheat: Error; no water vapor tracer !!!"
196            write(*,*) "H2O is always needed if calleuv=.true."
197            stop
198         else
199            nespeuv=nespeuv+1
200         endif
201         g_o1d=igcm_o1d
202         if (g_o1d.eq.0) then
203            write(*,*) "euvheat: Error; no O1D tracer !!!"
204            write(*,*) "O1D must always be present if thermochem=T"
205            stop
206         else
207            nespeuv=nespeuv+1 
208         endif
209         g_co=igcm_co
210         if (g_co.eq.0) then
211            write(*,*) "euvheat: Error; no CO tracer !!!"
212            write(*,*) "CO is always needed if calleuv=.true."
213            stop
214         else
215            nespeuv=nespeuv+1
216         endif
217         g_h=igcm_h
218         if (g_h.eq.0) then
219            write(*,*) "euvheat: Error; no H tracer !!!"
220            write(*,*) "H is always needed if calleuv=.true."
221            stop
222         else
223            nespeuv=nespeuv+1
224         endif
225         
226         euvmod = 0             !Default: C/O/H chemistry
227         !Check if O3 is present
228         g_o3=igcm_o3
229         if (g_o3.eq.0) then
230            write(*,*) "euvheat: Error; no O3 tracer !!!"
231            write(*,*) "O3 must be present if calleuv=.true."
232            stop
233         else
234            nespeuv=nespeuv+1
235            euvmod=1
236         endif
237
238         !Nitrogen species
239         !NO is used to determine if N chemistry is wanted
240         !euvmod=2 -> N chemistry
241         g_no=igcm_no
242         if (g_no.eq.0) then
243            write(*,*) "euvheat: no NO tracer"
244            write(*,*) "No N species in UV heating"
245         else if(g_no.ne.0) then
246            nespeuv=nespeuv+1
247            euvmod=2
248         endif
249         ! n2
250         g_n2=igcm_n2
251         if(euvmod.eq.2) then
252            if (g_n2.eq.0) then
253               write(*,*) "euvheat: Error; no N2 tracer !!!"
254               write(*,*) "N2 needed if NO is in traceur.def"
255               stop
256            else
257               nespeuv=nespeuv+1
258            endif
259         endif  ! Of if(euvmod.eq.2)
260         ! N
261         g_n=igcm_n
262         if(euvmod == 2) then
263            if (g_n.eq.0) then
264               write(*,*) "euvheat: Error; no N tracer !!!"
265               write(*,*) "N needed if NO is in traceur.def"
266               stop
267            else if(g_n.ne.0) then
268               nespeuv=nespeuv+1
269            endif
270         else
271            if(g_n /= 0) then
272               write(*,*) "euvheat: Error: N present, but NO not!!!"
273               write(*,*) "Both must be in traceur.def"
274               stop
275            endif
276         endif   !Of if(euvmod==2)
277         !NO2
278         g_no2=igcm_no2
279         if(euvmod == 2) then
280            if (g_no2.eq.0) then
281               write(*,*) "euvheat: Error; no NO2 tracer !!!"
282               write(*,*) "NO2 needed if NO is in traceur.def"
283               stop
284            else if(g_no2.ne.0) then
285               nespeuv=nespeuv+1
286            endif
287         else
288            if(g_no2 /= 0) then
289               write(*,*) "euvheat: Error: NO2 present, but NO not!!!"
290               write(*,*) "Both must be in traceur.def"
291               stop
292            endif
293         endif   !Of if(euvmod==2)
294         !N2D
295         g_n2d=igcm_n2d
296         if(euvmod == 2) then
297            if (g_n2d.eq.0) then
298               write(*,*) "euvheat: Error; no N2D tracer !!!"
299               write(*,*) "N2D needed if NO is in traceur.def"
300               stop
301            else
302               nespeuv=nespeuv+1 
303            endif
304         else
305            if(g_n2d /= 0) then
306               write(*,*) "euvheat: Error: N2D present, but NO not!!!"
307               write(*,*) "Both must be in traceur.def"
308               stop
309            endif
310         endif   !Of if(euvmod==2)
311
312         !Check if nespeuv is appropriate for the value of euvmod
313         select case(euvmod)
314         case(0)
315            if(nespeuv.ne.11) then
316               write(*,*)'euvheat: Wrong number of tracers!'
317               stop
318            else
319               write(*,*)'euvheat: Computing absorption by',nespeuv, &
320                    ' species'
321            endif
322         case(1)
323            if(nespeuv.ne.12) then
324               write(*,*)'euvheat: Wrong number of tracers!',nespeuv
325               stop
326            else
327               write(*,*)'euvheat: Computing absorption by',nespeuv,  &
328                    ' species'
329            endif
330         case(2)
331            if(nespeuv.ne.17) then
332               write(*,*)'euvheat: Wrong number of tracers!'
333               stop
334            else
335               write(*,*)'euvheat: Computing absorption by',nespeuv,  &
336                    ' species'
337            endif
338         end select
339       
340         !Allocate density vector
341         allocate(rm(nlayer,nespeuv))
342
343         firstcall= .false.
344      endif                     ! of if (firstcall)
345
346!cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
347
348
349      ! build local updated values of tracers and temperature
350      do l=1,nlayer
351        do ig=1,ngrid
352          ! chemical species
353          zq(ig,l,g_co2)=pq(ig,l,g_co2)+pdq(ig,l,g_co2)*ptimestep
354          zq(ig,l,g_o2)=pq(ig,l,g_o2)+pdq(ig,l,g_o2)*ptimestep
355          zq(ig,l,g_o)=pq(ig,l,g_o)+pdq(ig,l,g_o)*ptimestep
356          zq(ig,l,g_h2)=pq(ig,l,g_h2)+pdq(ig,l,g_h2)*ptimestep
357          zq(ig,l,g_h2o2)=pq(ig,l,g_h2o2)+pdq(ig,l,g_h2o2)*ptimestep
358          zq(ig,l,g_h2o)=pq(ig,l,g_h2o)+pdq(ig,l,g_h2o)*ptimestep
359          zq(ig,l,g_n2)=pq(ig,l,g_n2)+pdq(ig,l,g_n2)*ptimestep
360          zq(ig,l,g_co)=pq(ig,l,g_co)+pdq(ig,l,g_co)*ptimestep
361          zq(ig,l,g_h)=pq(ig,l,g_h)+pdq(ig,l,g_h)*ptimestep
362          if(euvmod.ge.1)   &
363               zq(ig,l,g_o3)=pq(ig,l,g_o3)+pdq(ig,l,g_o3)*ptimestep
364          if(euvmod.eq.2) then
365             zq(ig,l,g_n)=pq(ig,l,g_n)+pdq(ig,l,g_n)*ptimestep
366             zq(ig,l,g_no)=pq(ig,l,g_no)+pdq(ig,l,g_no)*ptimestep
367             zq(ig,l,g_no2)=pq(ig,l,g_no2)+pdq(ig,l,g_no2)*ptimestep
368          endif
369          if(euvmod.gt.2.or.euvmod.lt.0) then
370             write(*,*)'euvheat: bad value for euvmod. Stop'
371             stop
372          endif
373          ! atmospheric temperature
374          zt(ig,l)=pt(ig,l)+pdt(ig,l)*ptimestep
375        enddo
376      enddo
377     
378     
379      do ig=1,ngrid
380         zenit=acos(mu0(ig))*180./acos(-1.)
381         
382         do l=1,nlayer
383            !Conversion to number density
384            dens=pplay(ig,l)/(rnew(ig,l)*zt(ig,l)) / 1.66e-21
385            rm(l,i_co2)  = zq(ig,l,g_co2)  * dens / mmol(g_co2)
386            rm(l,i_o2)   = zq(ig,l,g_o2)   * dens / mmol(g_o2)
387            rm(l,i_o)    = zq(ig,l,g_o)    * dens / mmol(g_o)
388            rm(l,i_h2)   = zq(ig,l,g_h2)   * dens / mmol(g_h2)
389            rm(l,i_h2o)  = zq(ig,l,g_h2o)  * dens / mmol(g_h2o)
390            rm(l,i_h2o2) = zq(ig,l,g_h2o2) * dens / mmol(g_h2o2)
391            rm(l,i_co)   = zq(ig,l,g_co)   * dens / mmol(g_co)
392            rm(l,i_h)    = zq(ig,l,g_h)    * dens / mmol(g_h)
393            !Only if O3, N or ion chemistry requested
394            if(euvmod.ge.1)   &
395                 rm(l,i_o3)   = zq(ig,l,g_o3)   * dens / mmol(g_o3)
396            !Only if N or ion chemistry requested
397            if(euvmod.ge.2) then
398               rm(l,i_n2)   = zq(ig,l,g_n2)   * dens / mmol(g_n2)
399               rm(l,i_n)    = zq(ig,l,g_n)    * dens / mmol(g_n)
400               rm(l,i_no)   = zq(ig,l,g_no)   * dens / mmol(g_no)         
401               rm(l,i_no2)  = zq(ig,l,g_no2)  * dens / mmol(g_no2)
402            endif
403         enddo
404
405!        zlocal(1)=-log(pplay(ig,1)/pplev(ig,1))
406!     &            *Rnew(ig,1)*zt(ig,1)/g
407         zlocal(1)=zzlay(ig,1)
408         zlocal(1)=zlocal(1)/1000.
409         tx(1)=zt(ig,1)
410
411         do l=2,nlayer
412            tx(l)=zt(ig,l)
413            zlocal(l)=zzlay(ig,l)/1000.
414         enddo
415        !Routine to calculate the UV heating
416         call hrtherm (ig,nlayer,euvmod,rm,nespeuv,tx,zlocal,zenit,zday,jtot)
417
418!        euveff=0.16    !UV heating efficiency. Following Fox et al. ASR 1996
419                       !should vary between 19% and 23%. Lower values
420                       !(i.e. 16%) can be used to compensate underestimation
421                       !of 15-um cooling (see Forget et al. JGR 2009 and
422                       !Gonzalez-Galindo et al. JGR 2009) for details
423        !Calculates the UV heating from the total photoabsorption coefficient
424        do l=1,nlayer
425          pdteuv(ig,l)=euveff*jtot(l)/10.                  &
426               /(cpnew(ig,l)*pplay(ig,l)/(rnew(ig,l)*zt(ig,l)))
427!     &                  *(1.52/dist_sol)**2  !The solar flux calculated in
428                                              !flujo.F is already corrected for
429                                              !the actual Mars-Sun distance
430        enddo   
431      enddo  ! of do ig=1,ngrid
432       
433      !Deallocations
434!      deallocate(rm)
435
436      return
437      end
Note: See TracBrowser for help on using the repository browser.