source: trunk/LMDZ.VENUS/libf/phyvenus/euvheat.F90 @ 3884

Last change on this file since 3884 was 3884, checked in by ikovalenko, 4 months ago
File size: 14.3 KB
Line 
1      SUBROUTINE euvheat(nlon, nlev,nqmx, pt,pplev,pplay,zzlay, &
2                   mu0,ptimestep,ptime,pq, pdq, pdteuv)
3
4        use chemparam_mod
5        use dimphy
6        use conc, only:  rnew, cpnew
7        use mmol_mod
8        use clesphys_mod
9        use YOMCST_mod
10
11      IMPLICIT NONE
12!=======================================================================
13!   subject:
14!   --------
15!   Computing heating rate due to EUV absorption
16!
17!   author:  MAC 2002
18!   ------
19!
20!   input:
21!   -----
22!   mu0(klon)           
23!   pplay(ngrid,nlayer)   pressure at middle of layers (Pa)
24!   zzlay                 ! altitude at the middle of the layers (m)
25!
26!   output:
27!   -------
28!
29!   pdteuv(ngrid,nlayer)      Heating rate (K/s)
30!
31!=======================================================================
32!
33!    0.  Declarations :
34!    ------------------
35!
36!#include "YOMCST.h"
37!#include "clesphys.h"
38!#include "mmol.h"
39!-----------------------------------------------------------------------
40!    Input/Output
41!    ------------
42
43
44      integer :: nlon
45      integer :: nlev
46      integer :: nqmx
47
48      real :: pt(nlon,nlev)
49      real :: pplev(nlon,nlev+1)
50      real :: pplay(nlon,nlev)
51      real :: zzlay(nlon,nlev)
52      real :: mu0(nlon)
53
54      real :: ptimestep,ptime
55      real :: pq(nlon,nlev,nqmx)
56      real :: pdq(nlon,nlev,nqmx)
57      real :: pdteuv(nlon,nlev)
58!
59!    Local variables :
60!    -----------------
61
62      integer,save :: nespeuv=17    ! Number of species considered (11, 12 or 17 (with nitrogen))
63      integer,save :: nspeuv_vgcm    ! Number of species considered currently considered into VGCM
64
65
66      INTEGER :: l,ig,n
67      integer,save :: euvmod = 0     !0: 4 (main) species  1: O3 chemistry 2: N chemistry, 3: C/O/H
68      real, allocatable, save :: rm(:,:)   !  number density (cm-3)
69      real :: zq(nlon,nlev,nqmx) ! local updated tracer quantity
70      real :: zt(nlon,nlev)      ! local updated atmospheric temperature
71      real :: zlocal(nlev)
72      real :: zenit
73      real :: jtot(nlev)
74      real :: dens              ! Total number density (cm-3)
75      real :: tx(nlev)
76       
77! tracer indexes for the EUV heating:
78!!! ATTENTION. These values have to be identical to those in hrterm.F90
79!!! If the values are changed there, the same has to be done here  !!!
80     
81!      integer,parameter :: ix_co2=1
82!      integer,parameter :: ix_o=3
83!      integer,parameter :: ix_co=4
84!      integer,parameter :: ix_n2=13
85
86      integer,parameter :: ix_co2  =  1
87      integer,parameter :: ix_co   =  2
88      integer,parameter :: ix_o    =  3
89      integer,parameter :: ix_o1d  =  4
90      integer,parameter :: ix_o2   =  5
91      integer,parameter :: ix_o3   =  6
92      integer,parameter :: ix_h    =  7
93      integer,parameter :: ix_h2   =  8
94      integer,parameter :: ix_oh   =  9
95      integer,parameter :: ix_ho2  = 10
96      integer,parameter :: ix_h2o2 = 11
97      integer,parameter :: ix_h2o  = 12
98      integer,parameter :: ix_n    = 13
99      integer,parameter :: ix_n2d  = 14
100      integer,parameter :: ix_no   = 15
101      integer,parameter :: ix_no2  = 16
102      integer,parameter :: ix_n2   = 17
103
104! Tracer indexes in the GCM:
105      integer,save :: g_co2=0
106      integer,save :: g_o=0
107      integer,save :: g_o2=0
108      integer,save :: g_h2=0
109      integer,save :: g_h2o2=0
110      integer,save :: g_h2o=0
111      integer,save :: g_o3=0
112      integer,save :: g_n2=0
113      integer,save :: g_n=0
114      integer,save :: g_no=0
115      integer,save :: g_co=0
116      integer,save :: g_h=0
117      integer,save :: g_no2=0
118      integer,save :: g_oh=0
119      integer,save :: g_ho2=0
120      integer,save :: g_o1d=0
121      integer,save :: g_n2d=0
122     
123      logical,save :: firstcall=.true.
124
125! Initializations and sanity checks:
126
127
128      if (firstcall) then
129         nspeuv_vgcm=0
130!        ! identify the indexes of the tracers we'll need
131         g_co2=i_co2
132         if (g_co2.eq.0) then
133            write(*,*) "euvheat: Error; no CO2 tracer !!!"
134            write(*,*) "CO2 is always needed if calleuv=.true."
135            stop
136         else
137            nspeuv_vgcm=nspeuv_vgcm+1
138         endif
139         g_o=i_o
140         if (g_o.eq.0) then
141            write(*,*) "euvheat: Error; no O tracer !!!"
142!            write(*,*) "O is always needed if calleuv=.true."
143            stop
144         else
145            nspeuv_vgcm=nspeuv_vgcm+1
146         endif
147         g_co=i_co
148         if (g_co.eq.0) then
149            write(*,*) "euvheat: Error; no CO tracer !!!"
150!            write(*,*) "CO is always needed if calleuv=.true."
151            stop
152         else
153            nspeuv_vgcm=nspeuv_vgcm+1
154         endif
155         ! n2
156         g_n2=i_n2
157         if (g_n2.eq.0) then
158            write(*,*) "euvheat: Error; no N2 tracer !!!"
159!                write(*,*) "N2 needed if NO is in traceur.def"
160            stop
161         else
162            nspeuv_vgcm=nspeuv_vgcm+1
163         endif
164         g_o2=i_o2
165         if (g_o2.eq.0) then
166            write(*,*) "euvheat: Error; no O2 tracer !!!"
167!            write(*,*) "O2 is always needed if calleuv=.true."
168            stop
169         else
170            nspeuv_vgcm=nspeuv_vgcm+1
171         endif
172         g_h2=i_h2
173         if (g_h2.eq.0) then
174            write(*,*) "euvheat: Error; no H2 tracer !!!"
175!            write(*,*) "H2 is always needed if calleuv=.true."
176            stop
177         else
178            nspeuv_vgcm=nspeuv_vgcm+1
179         endif
180         g_oh=i_oh
181         if (g_oh.eq.0) then
182            write(*,*) "euvheat: Error; no OH tracer !!!"
183!            write(*,*) "OH must always be present if thermochem=T"
184            stop
185         else
186            nspeuv_vgcm=nspeuv_vgcm+1 
187         endif
188         g_ho2=i_ho2
189         if (g_ho2.eq.0) then
190            write(*,*) "euvheat: Error; no HO2 tracer !!!"
191!            write(*,*) "HO2 must always be present if thermochem=T"
192            stop
193         else
194            nspeuv_vgcm=nspeuv_vgcm+1 
195         endif
196         g_h2o2=i_h2o2
197         if (g_h2o2.eq.0) then
198            write(*,*) "euvheat: Error; no H2O2 tracer !!!"
199!            write(*,*) "H2O2 is always needed if calleuv=.true."
200            stop
201         else
202            nspeuv_vgcm=nspeuv_vgcm+1
203         endif
204         g_h2o=i_h2o
205         if (g_h2o.eq.0) then
206            write(*,*) "euvheat: Error; no water vapor tracer !!!"
207!            write(*,*) "H2O is always needed if calleuv=.true."
208            stop
209         else
210            nspeuv_vgcm=nspeuv_vgcm+1
211         endif
212         g_o1d=i_o1d
213         if (g_o1d.eq.0) then
214            write(*,*) "euvheat: Error; no O1D tracer !!!"
215!            write(*,*) "O1D must always be present if thermochem=T"
216            stop
217         else
218            nspeuv_vgcm=nspeuv_vgcm+1 
219         endif
220         g_h=i_h
221         if (g_h.eq.0) then
222            write(*,*) "euvheat: Error; no H tracer !!!"
223!            write(*,*) "H is always needed if calleuv=.true."
224            stop
225         else
226            nspeuv_vgcm=nspeuv_vgcm+1
227         endif
228         
229!         euvmod = 1            !Default: C/O/H chemistry
230!         !Check if O3 is present
231         g_o3=i_o3
232         if (g_o3.eq.0) then
233            write(*,*) "euvheat: Error; no O3 tracer !!!"
234            write(*,*) "O3 must be present if calleuv=.true."
235            stop
236         else
237            nspeuv_vgcm=nspeuv_vgcm+1
238            euvmod=1
239         endif
240
241         !Nitrogen species
242         !NO is used to determine if N chemistry is wanted
243         !euvmod=2 -> N chemistry
244         g_no=i_no
245         if (g_no.eq.0) then
246            write(*,*) "euvheat: no NO tracer"
247            write(*,*) "No N species in UV heating"
248         else if(g_no.ne.0) then
249            nspeuv_vgcm=nspeuv_vgcm+1
250            euvmod=2
251         endif
252         ! N
253         g_n=i_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               nspeuv_vgcm=nspeuv_vgcm+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=i_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               nspeuv_vgcm=nspeuv_vgcm+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=i_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               nspeuv_vgcm=nspeuv_vgcm+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
333         !Allocate density vector
334         allocate(rm(nlev,nespeuv))
335
336         firstcall= .false.
337      endif                     ! of if (firstcall)
338
339!      write(*,*),  "CHECK n species currently used into VGCM",  nspeuv_vgcm
340
341
342!cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
343
344      ! euvmod selection security
345      if(euvmod.gt.2.or.euvmod.lt.0) then
346         write(*,*)'euvheat: bad value for euvmod. Stop'
347         stop
348      endif
349
350      ! build local updated values of tracers (if any) and temperature
351      do l=1,nlev
352        do ig=1,nlon
353          ! chemical species
354          zq(ig,l,g_co2)=pq(ig,l,g_co2)     ! CO2 
355          zq(ig,l,g_o2)=pq(ig,l,g_o2)       ! O2
356          zq(ig,l,g_o)=pq(ig,l,g_o)         ! O(3P)
357          zq(ig,l,g_h2)=pq(ig,l,g_h2)       ! H2
358          zq(ig,l,g_h2o2)=pq(ig,l,g_h2o2)   ! H2O2
359          zq(ig,l,g_h2o)=pq(ig,l,g_h2o)     ! H2O
360          zq(ig,l,g_n2)=pq(ig,l,g_n2)       ! N2   
361          zq(ig,l,g_co)=pq(ig,l,g_co)       ! CO
362          zq(ig,l,g_h)=pq(ig,l,g_h)         ! H
363
364          !Only if O3, N or ion chemistry requested
365          if(euvmod.ge.1) then
366             zq(ig,l,g_o3)=pq(ig,l,g_o3)    ! 03
367          endif
368          !Only if N or ion chemistry requested
369          if(euvmod.ge.2) then
370             zq(ig,l,g_n)=pq(ig,l,g_n)      ! N
371             zq(ig,l,g_no)=pq(ig,l,g_no)    ! NO
372             zq(ig,l,g_no2)=pq(ig,l,g_no2)  ! NO2
373          endif
374          ! atmospheric temperature
375          zt(ig,l)=pt(ig,l)
376!      write(*,*),  "CHECK update densities L332 euv",   zq(ig,l,g_co2)
377        enddo
378      enddo
379     
380      !Solar flux calculation     
381      do ig=1,nlon
382         zenit=acos(mu0(ig))*180./acos(-1.)  !convers from rad to deg
383         
384         do l=1,nlev         
385            !Conversion to number density   
386            !!!  use R specific = R/MolarMass
387            dens=pplay(ig,l)/(rnew(ig,l)*zt(ig,l)) / 1.66e-21   ! [g mol-1] [cm-3]         
388
389            rm(l,ix_co2)  = zq(ig,l,g_co2) * dens / M_tr(g_co2) ! [cm-3]
390            rm(l,ix_o2)   = zq(ig,l,g_o2)  * dens / M_tr(g_o2)
391            rm(l,ix_o)    = zq(ig,l,g_o)   * dens / M_tr(g_o)
392            rm(l,ix_h2)   = zq(ig,l,g_h2)  * dens / M_tr(g_h2)
393            rm(l,ix_h2o)  = zq(ig,l,g_h2o) * dens / M_tr(g_h2o)
394            rm(l,ix_h2o2) = zq(ig,l,g_h2o2)* dens / M_tr(g_h2o2)
395            rm(l,ix_co)   = zq(ig,l,g_co)  * dens / M_tr(g_co)
396            rm(l,ix_n2)   = zq(ig,l,g_n2)  * dens / M_tr(g_n2)
397            rm(l,ix_h)    = zq(ig,l,g_h)   * dens / M_tr(g_h)
398
399            !Only if O3, N or ion chemistry requested
400            if(euvmod.ge.1) then
401               rm(l,ix_o3)   = zq(ig,l,g_o3) * dens / M_tr(g_o3)
402            endif
403            !Only if N or ion chemistry requested
404            if(euvmod.ge.2) then
405               rm(l,ix_n)    = zq(ig,l,g_n)    * dens / M_tr(g_n)
406               rm(l,ix_no)   = zq(ig,l,g_no)   * dens / M_tr(g_no)
407               rm(l,ix_no2)  = zq(ig,l,g_no2)  * dens / M_tr(g_no2)
408            endif
409
410!           write(*,*),  "CHECK n density", l, rm(l,ix_co2)
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.    ! conversion m ---> km
417         tx(1)=zt(ig,1)
418         
419         do l=2,nlev
420            tx(l)=zt(ig,l)
421            zlocal(l)=zzlay(ig,l)/1000.
422         enddo
423
424        !Routine to calculate the UV heating
425         call hrtherm (ig,euvmod,rm,nespeuv,tx,zlocal,zenit,jtot)
426         
427
428                       !euveff: 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
431                       !underestimation
432                       !of 15-um cooling (see Forget et al. JGR 2009 and
433                       !Gonzalez-Galindo et al. JGR 2009) for details
434
435        !Calculates the UV heating from the total photoabsorption coefficient
436        do l=1,nlev
437!       jtot conversion from erg/(s*cm3) ---> J/(s*m3)
438          pdteuv(ig,l)=euveff*jtot(l)/10.                  &
439               /(cpnew(ig,l)*pplay(ig,l)/(rnew(ig,l)*zt(ig,l)))
440
441        enddo   
442      enddo  ! of do ig=1,nlon
443
444      !Deallocations
445      !deallocate(rm)
446
447      return
448      end
Note: See TracBrowser for help on using the repository browser.