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

Last change on this file since 3948 was 3726, checked in by emillour, 7 months ago

Mars PCM:
Turn "callkeys.h" into module "callkeys_mod.F90"
EM

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