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

Last change on this file since 1242 was 1226, checked in by aslmd, 11 years ago

LMDZ.MARS : Replaced comcstfi and planete includes by modules.

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