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

Last change on this file since 3567 was 3464, checked in by emillour, 2 months ago

Mars PCM:
Some tidying in aeronomars:

  • make a jthermcalc_util.F to contain utility routines (used by jthermcal & jthermcalc_e107). Also add the possibility (turned off by default) in the interfast routine to do extra sanity checks.
  • turn chemthermos, euvheat, hrtherm, jthermcalc, jthermcalc_e107, paramphoto_compact and photochemistry into modules.

EM

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