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

Last change on this file since 1119 was 1119, checked in by emillour, 11 years ago

Mars GCM:

  • Bug fix: Sun-Mars distance was not correctly taken into account in the solvarmod==1 (daily varying realistic EUV input) case.

FGG

File size: 14.6 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 "dimensions.h"
34!#include "dimphys.h"
35!#include "comcstfi.h"
36#include "callkeys.h"
37!#include "comdiurn.h"
38!#include "param.h"
39!#include "param_v4.h"
40!#include "chimiedata.h"
41!#include "tracer.h"
42!#include "conc.h"
43!-----------------------------------------------------------------------
44!    Input/Output
45!    ------------
46
47      integer,intent(in) :: ngrid ! number of atmospheric columns
48      integer,intent(in) :: nlayer ! number of atmospheric layers
49      integer,intent(in) :: nq ! number of advected tracers
50      real :: pt(ngrid,nlayer)
51      real :: pdt(ngrid,nlayer)
52      real :: pplev(ngrid,nlayer+1)
53      real :: pplay(ngrid,nlayer)
54      real :: zzlay(ngrid,nlayer)
55      real :: mu0(ngrid)
56      real :: ptimestep,ptime
57      real :: zday
58      real :: pq(ngrid,nlayer,nq)
59      real :: pdq(ngrid,nlayer,nq)
60
61      real :: pdteuv(ngrid,nlayer)
62!
63!    Local variables :
64!    -----------------
65      integer,save :: nespeuv    ! Number of species considered
66
67      INTEGER :: l,ig,n
68      integer,save :: euvmod
69      real, allocatable, save :: rm(:,:)   ! number density (cm-3)
70      real :: zq(ngrid,nlayer,nq) ! local updated tracer quantity
71      real :: zt(ngrid,nlayer)      ! local updated atmospheric temperature
72      real :: zlocal(nlayer)
73      real :: zenit
74      real :: jtot(nlayer)
75      real :: dens                                              ! amu/cm-3
76      real :: tx(nlayer)
77!      real euveff     !UV heating efficiency
78     
79! tracer indexes for the EUV heating:
80!!! ATTENTION. These values have to be identical to those in chemthermos.F90
81!!! If the values are changed there, the same has to be done here  !!!
82      integer,parameter :: i_co2=1
83      integer,parameter :: i_o2=2
84      integer,parameter :: i_o=3
85      integer,parameter :: i_co=4
86      integer,parameter :: i_h=5
87      integer,parameter :: i_oh=6
88      integer,parameter :: i_ho2=7
89      integer,parameter :: i_h2=8
90      integer,parameter :: i_h2o=9
91      integer,parameter :: i_h2o2=10
92      integer,parameter :: i_o1d=11
93      integer,parameter :: i_o3=12
94      integer,parameter :: i_n2=13
95      integer,parameter :: i_n=14
96      integer,parameter :: i_no=15
97      integer,parameter :: i_n2d=16
98      integer,parameter :: i_no2=17
99
100     
101! Tracer indexes in the GCM:
102      integer,save :: g_co2=0
103      integer,save :: g_o=0
104      integer,save :: g_o2=0
105      integer,save :: g_h2=0
106      integer,save :: g_h2o2=0
107      integer,save :: g_h2o=0
108      integer,save :: g_o3=0
109      integer,save :: g_n2=0
110      integer,save :: g_n=0
111      integer,save :: g_no=0
112      integer,save :: g_co=0
113      integer,save :: g_h=0
114      integer,save :: g_no2=0
115      integer,save :: g_oh=0
116      integer,save :: g_ho2=0
117      integer,save :: g_o1d=0
118      integer,save :: g_n2d=0
119
120
121      logical,save :: firstcall=.true.
122
123! Initializations and sanity checks:
124
125
126      if (firstcall) then
127         nespeuv=0
128        ! identify the indexes of the tracers we'll need
129         g_co2=igcm_co2
130         if (g_co2.eq.0) then
131            write(*,*) "euvheat: Error; no CO2 tracer !!!"
132            write(*,*) "CO2 is always needed if calleuv=.true."
133            stop
134         else
135            nespeuv=nespeuv+1
136         endif
137         g_o=igcm_o
138         if (g_o.eq.0) then
139            write(*,*) "euvheat: Error; no O tracer !!!"
140            write(*,*) "O is always needed if calleuv=.true."
141            stop
142         else
143            nespeuv=nespeuv+1
144         endif
145         g_o2=igcm_o2
146         if (g_o2.eq.0) then
147            write(*,*) "euvheat: Error; no O2 tracer !!!"
148            write(*,*) "O2 is always needed if calleuv=.true."
149            stop
150         else
151            nespeuv=nespeuv+1
152         endif
153         g_h2=igcm_h2
154         if (g_h2.eq.0) then
155            write(*,*) "euvheat: Error; no H2 tracer !!!"
156            write(*,*) "H2 is always needed if calleuv=.true."
157            stop
158         else
159            nespeuv=nespeuv+1
160         endif
161         g_oh=igcm_oh
162         if (g_oh.eq.0) then
163            write(*,*) "euvheat: Error; no OH tracer !!!"
164            write(*,*) "OH must always be present if thermochem=T"
165            stop
166         else
167            nespeuv=nespeuv+1 
168         endif
169         g_ho2=igcm_ho2
170         if (g_ho2.eq.0) then
171            write(*,*) "euvheat: Error; no HO2 tracer !!!"
172            write(*,*) "HO2 must always be present if thermochem=T"
173            stop
174         else
175            nespeuv=nespeuv+1 
176         endif
177         g_h2o2=igcm_h2o2
178         if (g_h2o2.eq.0) then
179            write(*,*) "euvheat: Error; no H2O2 tracer !!!"
180            write(*,*) "H2O2 is always needed if calleuv=.true."
181            stop
182         else
183            nespeuv=nespeuv+1
184         endif
185         g_h2o=igcm_h2o_vap
186         if (g_h2o.eq.0) then
187            write(*,*) "euvheat: Error; no water vapor tracer !!!"
188            write(*,*) "H2O is always needed if calleuv=.true."
189            stop
190         else
191            nespeuv=nespeuv+1
192         endif
193         g_o1d=igcm_o1d
194         if (g_o1d.eq.0) then
195            write(*,*) "euvheat: Error; no O1D tracer !!!"
196            write(*,*) "O1D must always be present if thermochem=T"
197            stop
198         else
199            nespeuv=nespeuv+1 
200         endif
201         g_co=igcm_co
202         if (g_co.eq.0) then
203            write(*,*) "euvheat: Error; no CO tracer !!!"
204            write(*,*) "CO is always needed if calleuv=.true."
205            stop
206         else
207            nespeuv=nespeuv+1
208         endif
209         g_h=igcm_h
210         if (g_h.eq.0) then
211            write(*,*) "euvheat: Error; no H tracer !!!"
212            write(*,*) "H is always needed if calleuv=.true."
213            stop
214         else
215            nespeuv=nespeuv+1
216         endif
217         
218         euvmod = 0             !Default: C/O/H chemistry
219         !Check if O3 is present
220         g_o3=igcm_o3
221         if (g_o3.eq.0) then
222            write(*,*) "euvheat: Error; no O3 tracer !!!"
223            write(*,*) "O3 must be present if calleuv=.true."
224            stop
225         else
226            nespeuv=nespeuv+1
227            euvmod=1
228         endif
229
230         !Nitrogen species
231         !NO is used to determine if N chemistry is wanted
232         !euvmod=2 -> N chemistry
233         g_no=igcm_no
234         if (g_no.eq.0) then
235            write(*,*) "euvheat: no NO tracer"
236            write(*,*) "No N species in UV heating"
237         else if(g_no.ne.0) then
238            nespeuv=nespeuv+1
239            euvmod=2
240         endif
241         ! n2
242         g_n2=igcm_n2
243         if(euvmod.eq.2) then
244            if (g_n2.eq.0) then
245               write(*,*) "euvheat: Error; no N2 tracer !!!"
246               write(*,*) "N2 needed if NO is in traceur.def"
247               stop
248            else
249               nespeuv=nespeuv+1
250            endif
251         endif  ! Of if(euvmod.eq.2)
252         ! N
253         g_n=igcm_n
254         if(euvmod == 2) then
255            if (g_n.eq.0) then
256               write(*,*) "euvheat: Error; no N tracer !!!"
257               write(*,*) "N needed if NO is in traceur.def"
258               stop
259            else if(g_n.ne.0) then
260               nespeuv=nespeuv+1
261            endif
262         else
263            if(g_n /= 0) then
264               write(*,*) "euvheat: Error: N present, but NO not!!!"
265               write(*,*) "Both must be in traceur.def"
266               stop
267            endif
268         endif   !Of if(euvmod==2)
269         !NO2
270         g_no2=igcm_no2
271         if(euvmod == 2) then
272            if (g_no2.eq.0) then
273               write(*,*) "euvheat: Error; no NO2 tracer !!!"
274               write(*,*) "NO2 needed if NO is in traceur.def"
275               stop
276            else if(g_no2.ne.0) then
277               nespeuv=nespeuv+1
278            endif
279         else
280            if(g_no2 /= 0) then
281               write(*,*) "euvheat: Error: NO2 present, but NO not!!!"
282               write(*,*) "Both must be in traceur.def"
283               stop
284            endif
285         endif   !Of if(euvmod==2)
286         !N2D
287         g_n2d=igcm_n2d
288         if(euvmod == 2) then
289            if (g_n2d.eq.0) then
290               write(*,*) "euvheat: Error; no N2D tracer !!!"
291               write(*,*) "N2D needed if NO is in traceur.def"
292               stop
293            else
294               nespeuv=nespeuv+1 
295            endif
296         else
297            if(g_n2d /= 0) then
298               write(*,*) "euvheat: Error: N2D present, but NO not!!!"
299               write(*,*) "Both must be in traceur.def"
300               stop
301            endif
302         endif   !Of if(euvmod==2)
303
304         !Check if nespeuv is appropriate for the value of euvmod
305         select case(euvmod)
306         case(0)
307            if(nespeuv.ne.11) then
308               write(*,*)'euvheat: Wrong number of tracers!'
309               stop
310            else
311               write(*,*)'euvheat: Computing absorption by',nespeuv, &
312                    ' species'
313            endif
314         case(1)
315            if(nespeuv.ne.12) then
316               write(*,*)'euvheat: Wrong number of tracers!',nespeuv
317               stop
318            else
319               write(*,*)'euvheat: Computing absorption by',nespeuv,  &
320                    ' species'
321            endif
322         case(2)
323            if(nespeuv.ne.17) then
324               write(*,*)'euvheat: Wrong number of tracers!'
325               stop
326            else
327               write(*,*)'euvheat: Computing absorption by',nespeuv,  &
328                    ' species'
329            endif
330         end select
331       
332         !Allocate density vector
333         allocate(rm(nlayer,nespeuv))
334
335         firstcall= .false.
336      endif                     ! of if (firstcall)
337
338!cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
339
340
341      ! build local updated values of tracers and temperature
342      do l=1,nlayer
343        do ig=1,ngrid
344          ! chemical species
345          zq(ig,l,g_co2)=pq(ig,l,g_co2)+pdq(ig,l,g_co2)*ptimestep
346          zq(ig,l,g_o2)=pq(ig,l,g_o2)+pdq(ig,l,g_o2)*ptimestep
347          zq(ig,l,g_o)=pq(ig,l,g_o)+pdq(ig,l,g_o)*ptimestep
348          zq(ig,l,g_h2)=pq(ig,l,g_h2)+pdq(ig,l,g_h2)*ptimestep
349          zq(ig,l,g_h2o2)=pq(ig,l,g_h2o2)+pdq(ig,l,g_h2o2)*ptimestep
350          zq(ig,l,g_h2o)=pq(ig,l,g_h2o)+pdq(ig,l,g_h2o)*ptimestep
351          zq(ig,l,g_n2)=pq(ig,l,g_n2)+pdq(ig,l,g_n2)*ptimestep
352          zq(ig,l,g_co)=pq(ig,l,g_co)+pdq(ig,l,g_co)*ptimestep
353          zq(ig,l,g_h)=pq(ig,l,g_h)+pdq(ig,l,g_h)*ptimestep
354          if(euvmod.ge.1)   &
355               zq(ig,l,g_o3)=pq(ig,l,g_o3)+pdq(ig,l,g_o3)*ptimestep
356          if(euvmod.eq.2) then
357             zq(ig,l,g_n)=pq(ig,l,g_n)+pdq(ig,l,g_n)*ptimestep
358             zq(ig,l,g_no)=pq(ig,l,g_no)+pdq(ig,l,g_no)*ptimestep
359             zq(ig,l,g_no2)=pq(ig,l,g_no2)+pdq(ig,l,g_no2)*ptimestep
360          endif
361          if(euvmod.gt.2.or.euvmod.lt.0) then
362             write(*,*)'euvheat: bad value for euvmod. Stop'
363             stop
364          endif
365          ! atmospheric temperature
366          zt(ig,l)=pt(ig,l)+pdt(ig,l)*ptimestep
367        enddo
368      enddo
369     
370      !Solar flux calculation
371      call flujo(solarcondate)
372
373                                         ! Not recommended for long runs
374                                         ! (e.g. to build the MCD), as the
375                                         ! solar conditions at the end will
376                                         ! be different to the ones initially
377                                         ! set
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,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.