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

Last change on this file since 3026 was 2615, checked in by romain.vande, 3 years ago

LMDZ_MARS RV : Open_MP;
Put all the "save" variables as "!$OMP THREADPRIVATE" in aeronomars

File size: 14.9 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)=-log(pplay(ig,1)/pplev(ig,1))
414!     &            *Rnew(ig,1)*zt(ig,1)/g
415         zlocal(1)=zzlay(ig,1)
416         zlocal(1)=zlocal(1)/1000.
417         tx(1)=zt(ig,1)
418
419         do l=2,nlayer
420            tx(l)=zt(ig,l)
421            zlocal(l)=zzlay(ig,l)/1000.
422         enddo
423        !Routine to calculate the UV heating
424         call hrtherm (ig,nlayer,euvmod,rm,nespeuv,tx,zlocal,zenit,zday,jtot)
425
426
427
428!        euveff=0.16    !UV heating efficiency. Following Fox et al. ASR 1996
429                       !should vary between 19% and 23%. Lower values
430                       !(i.e. 16%) can be used to compensate underestimation
431                       !of 15-um cooling (see Forget et al. JGR 2009 and
432                       !Gonzalez-Galindo et al. JGR 2009) for details
433        !Calculates the UV heating from the total photoabsorption coefficient
434        do l=1,nlayer
435          pdteuv(ig,l)=euveff*jtot(l)/10.                  &
436               /(cpnew(ig,l)*pplay(ig,l)/(rnew(ig,l)*zt(ig,l)))
437!     &                  *(1.52/dist_sol)**2  !The solar flux calculated in
438                                              !flujo.F is already corrected for
439                                              !the actual Mars-Sun distance
440        enddo   
441      enddo  ! of do ig=1,ngrid
442       
443      !Deallocations
444!      deallocate(rm)
445
446      return
447      end
Note: See TracBrowser for help on using the repository browser.