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

Last change on this file since 3289 was 3158, checked in by csegonne, 13 months ago

MARS PCM
Cleaning of conduction.F, euvheat.F90, moldiff.F and molvis.F, some commented lines referring to a local calculation of layers/levels altitudes have been removed.

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