source: trunk/LMDZ.MARS/libf/aeronomars/chemthermos.F90

Last change on this file was 3464, checked in by emillour, 2 months ago

Mars PCM:
Some tidying in aeronomars:

  • make a jthermcalc_util.F to contain utility routines (used by jthermcal & jthermcalc_e107). Also add the possibility (turned off by default) in the interfast routine to do extra sanity checks.
  • turn chemthermos, euvheat, hrtherm, jthermcalc, jthermcalc_e107, paramphoto_compact and photochemistry into modules.

EM

File size: 20.6 KB
RevLine 
[3464]1  MODULE chemthermos_mod
2 
3  IMPLICIT NONE
4 
5  CONTAINS
6 
[1266]7      SUBROUTINE chemthermos(ig,nlayer,lswitch,chemthermod,zycol,ztemp, &
[1888]8           zdens,zpress,zlocal,zenit,ptimestep,zday,em_no,em_o2)
[635]9
[1036]10      use tracer_mod, only: nqmx, 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,  &
13                            igcm_n, igcm_no, igcm_no2, igcm_n2d,        &
14                            igcm_co2plus, igcm_o2plus, igcm_coplus,     &
15                            igcm_cplus, igcm_nplus, igcm_noplus,        &
16                            igcm_n2plus, igcm_hplus, igcm_hco2plus,     &
17                            igcm_elec, igcm_oplus
[1888]18
19      use param_v4_h, only: Pno, Po2
[3464]20      use jthermcalc_e107_mod, only: jthermcalc_e107
21      use paramfoto_compact_mod, only: paramfoto_compact
22
[635]23      IMPLICIT NONE
24!=======================================================================
25!   subject:
26!   --------
27!   Computing chemical variations in the thermosphere
28!
29!   author:  MAC July 2003
30!   ------
31!   modifications:
32!   -------------
33!     Ehouarn Sept 2008: added handling of tracers by their names
34!     Francisco Sept 2009: added ionosphere and N chemistry
35!     Francisco Mar 2012: Different chemical schemes possible according to
36!                         traceur.def
37!=======================================================================
38!
39!    0.  Declarations :
40!    ------------------
41!
[3464]42      include "callkeys.h"
[635]43!-----------------------------------------------------------------------
44!    Input/Output
45!    ------------
[1266]46      integer :: lswitch,ig,chemthermod,nlayer
47      real :: zycol(nlayer,nqmx)
48      real :: ztemp(nlayer)
49      real :: zdens(nlayer)
50      real :: zpress(nlayer)                    ! in mbar
51      real :: zlocal(nlayer)
[635]52      real :: zenit
53      real :: ptimestep
54      real :: zday
[1888]55      real :: em_no(nlayer)
56      real :: em_o2(nlayer)
[635]57!
58!    Local variables :
59!    -----------------
[705]60      integer :: l
[635]61      integer,save :: nesptherm
62      real,allocatable :: rm(:,:)               !number density (cm-3)
63      logical,save :: firstcall=.true.
64
[2615]65!$OMP THREADPRIVATE(nesptherm,firstcall)
66
[635]67! Tracer indexes in the GCM:
68      integer,save :: g_co2=0
69      integer,save :: g_co=0
70      integer,save :: g_o=0
71      integer,save :: g_o1d=0
72      integer,save :: g_o2=0
73      integer,save :: g_o3=0
74      integer,save :: g_h=0
75      integer,save :: g_h2=0
76      integer,save :: g_oh=0
77      integer,save :: g_ho2=0
78      integer,save :: g_h2o2=0
79      integer,save :: g_n2=0
80      integer,save :: g_h2o_vap=0
81      integer,save :: g_n=0
82      integer,save :: g_no=0
83      integer,save :: g_no2=0
84      integer,save :: g_n2d=0
85      integer,save :: g_co2plus=0
86      integer,save :: g_oplus=0
87      integer,save :: g_o2plus=0
88      integer,save :: g_coplus=0
89      integer,save :: g_cplus=0
90      integer,save :: g_nplus=0
91      integer,save :: g_noplus=0
92      integer,save :: g_n2plus=0
93      integer,save :: g_hplus=0
94      integer,save :: g_hco2plus=0
95      integer,save :: g_elec=0
[2615]96
97!$OMP THREADPRIVATE(g_co2,g_co,g_o,g_o1d,g_o2,g_o3,g_h,g_h2,g_oh,g_ho2,g_h2o2,g_n2,g_h2o_vap)
98!$OMP THREADPRIVATE(g_n,g_no,g_no2,g_n2d,g_co2plus,g_oplus,g_o2plus,g_coplus,g_cplus,g_nplus,g_noplus,g_n2plus,g_hplus,g_hco2plus,g_elec)
[635]99     
100! Tracer indexes in the thermospheric chemistry:
101!!! IMPORTANT. These values have to be the same than in
102!!! jthermcal.F, column.F, paramfoto_compact.F, euvheat.F
103!!! and hrtherm.F
104!!! If some value is changed here, the same change has to
105!!! be made into the other subroutines
[2042]106      integer,parameter :: i_co2  =  1
107      integer,parameter :: i_co   =  2
108      integer,parameter :: i_o    =  3
109      integer,parameter :: i_o1d  =  4
110      integer,parameter :: i_o2   =  5
111      integer,parameter :: i_o3   =  6
112      integer,parameter :: i_h    =  7
113      integer,parameter :: i_h2   =  8
114      integer,parameter :: i_oh   =  9
115      integer,parameter :: i_ho2  = 10
116      integer,parameter :: i_h2o2 = 11
117      integer,parameter :: i_h2o  = 12
118      integer,parameter :: i_n    = 13
119      integer,parameter :: i_n2d  = 14
120      integer,parameter :: i_no   = 15
121      integer,parameter :: i_no2  = 16
122      integer,parameter :: i_n2   = 17
123!      integer,parameter :: i_co2=1
124!      integer,parameter :: i_o2=2
125!      integer,parameter :: i_o=3
126!      integer,parameter :: i_co=4
127!      integer,parameter :: i_h=5
128!      integer,parameter :: i_oh=6
129!      integer,parameter :: i_ho2=7
130!      integer,parameter :: i_h2=8
131!      integer,parameter :: i_h2o=9
132!      integer,parameter :: i_h2o2=10
133!      integer,parameter :: i_o1d=11
134!      integer,parameter :: i_o3=12
135!      integer,parameter :: i_n2=13
136!      integer,parameter :: i_n=14
137!      integer,parameter :: i_no=15
138!      integer,parameter :: i_n2d=16
139!      integer,parameter :: i_no2=17
[635]140      integer,parameter :: i_co2plus=18
141      integer,parameter :: i_oplus=19
142      integer,parameter :: i_o2plus=20
143      integer,parameter :: i_coplus=21
144      integer,parameter :: i_cplus=22
145      integer,parameter :: i_nplus=23
146      integer,parameter :: i_noplus=24
147      integer,parameter :: i_n2plus=25
148      integer,parameter :: i_hplus=26
149      integer,parameter :: i_hco2plus=27
150      integer,parameter :: i_elec=28
151
152     
153! Initializations at first call
154      if (firstcall) then
155         !Number of species to consider
156         nesptherm=0
157         ! get the indexes of the tracers we'll need
158         ! It is not really necessary to check if the species are present,
159         ! as it is already done in calchim.F90. But we include it to double-check
160
161         ! Species always included if photochemistry
162         g_co2=igcm_co2
163         if (g_co2.eq.0) then
164            write(*,*) "chemthermos: Error; no CO2 tracer !!!"
165            write(*,*) "CO2 must always be present if thermochem=T"
166            stop
167         else
168            nesptherm=nesptherm+1 
169         endif
170         g_co=igcm_co
171         if (g_co.eq.0) then
172            write(*,*) "chemthermos: Error; no CO tracer !!!"
173            write(*,*) "CO must always be present if thermochem=T"
174            stop
175         else
176            nesptherm=nesptherm+1 
177         endif
178         g_o=igcm_o
179         if (g_o.eq.0) then
180            write(*,*) "chemthermos: Error; no O tracer !!!"
181            write(*,*) "O must always be present if thermochem=T"
182            stop
183         else
184            nesptherm=nesptherm+1 
185         endif
186         g_o1d=igcm_o1d
187         if (g_o1d.eq.0) then
188            write(*,*) "chemthermos: Error; no O1D tracer !!!"
189            write(*,*) "O1D must always be present if thermochem=T"
190            stop
191         else
192            nesptherm=nesptherm+1 
193         endif
194         g_o2=igcm_o2
195         if (g_o2.eq.0) then
196            write(*,*) "chemthermos: Error; no O2 tracer !!!"
197            write(*,*) "O2 must always be present if thermochem=T"
198            stop
199         else
200            nesptherm=nesptherm+1 
201         endif
202         g_h=igcm_h
203         if (g_h.eq.0) then
204            write(*,*) "chemthermos: Error; no H tracer !!!"
205            write(*,*) "H must always be present if thermochem=T"
206            stop
207         else
208            nesptherm=nesptherm+1 
209         endif
210         g_h2=igcm_h2
211         if (g_h2.eq.0) then
212            write(*,*) "chemthermos: Error; no H2 tracer !!!"
213            write(*,*) "H2 must always be present if thermochem=T"
214            stop
215         else
216            nesptherm=nesptherm+1 
217         endif
218         g_oh=igcm_oh
219         if (g_oh.eq.0) then
220            write(*,*) "chemthermos: Error; no OH tracer !!!"
221            write(*,*) "OH must always be present if thermochem=T"
222            stop
223         else
224            nesptherm=nesptherm+1 
225         endif
226         g_ho2=igcm_ho2
227         if (g_ho2.eq.0) then
228            write(*,*) "chemthermos: Error; no HO2 tracer !!!"
229            write(*,*) "HO2 must always be present if thermochem=T"
230            stop
231         else
232            nesptherm=nesptherm+1 
233         endif
234         g_h2o2=igcm_h2o2
235         if (g_h2o2.eq.0) then
236            write(*,*) "chemthermos: Error; no H2O2 tracer !!!"
237            write(*,*) "H2O2 must always be present if thermochem=T"
238            stop
239         else
240            nesptherm=nesptherm+1 
241         endif
242         g_h2o_vap=igcm_h2o_vap
243         if (g_h2o_vap.eq.0) then
244            write(*,*) "chemthermos: Error; no water vapor tracer !!!"
245            write(*,*) "H2O must always be present if thermochem=T"
246            stop
247         else
248            nesptherm=nesptherm+1 
249         endif
250         !Verify if O3 chemistry wanted
251         g_o3=igcm_o3
252         if(chemthermod.ge.1) then
253            if (g_o3.eq.0) then
254               write(*,*) "chemthermos: Error; no O3 tracer !!!"
255               write(*,*) "O3 must be present if NO is in traceur.def"
256               stop
257            else
258               nesptherm=nesptherm+1
259            endif
260         else
261            if(g_o3 == 0) then
262               write(*,*) "chemthermos: No O3 chemistry"
263            else               
264               write(*,*) "chemthermos: O3 chemistry included"
265               nesptherm=nesptherm+1 
266               chemthermod=1
267            endif
268         endif    ! Of if(chemthermod.ge.1)
269           
270         ! N chemistry
271         if(chemthermod.ge.2) then
272            g_n2=igcm_n2
273            if (g_n2.eq.0) then
274               write(*,*) "chemthermos: Error; no N2 tracer !!!"
275               write(*,*) "N2 must be present if NO is in traceur.def"
276               stop
277            else
278               nesptherm=nesptherm+1 
279            endif
280            g_n=igcm_n
281            if (g_n.eq.0) then
282               write(*,*) "chemthermos: Error; no N tracer !!!"
283               write(*,*) "N must be present if NO is in traceur.def"
284               stop
285            else
286               nesptherm=nesptherm+1 
287            endif
288            g_no=igcm_no
289            if (g_no.eq.0) then
290               write(*,*) "chemthermos: Error; no NO tracer !!!"
291               write(*,*) "NO must be present if N chemistry wanted"
292               stop
293            else
294               nesptherm=nesptherm+1 
295            endif
296            g_no2=igcm_no2
297            if (g_no2.eq.0) then
298               write(*,*) "chemthermos: Error; no NO2 tracer !!!"
299               write(*,*) "NO must be present if NO is in traceur.def"
300               stop
301            else
302               nesptherm=nesptherm+1 
303            endif
304            g_n2d=igcm_n2d
305            if (g_n2d.eq.0) then
306               write(*,*) "chemthermos: Error; no N2D tracer !!!"
307               write(*,*) "NO must be present if NO is in traceur.def"
308               stop
309            else
310               nesptherm=nesptherm+1 
311            endif
312         endif   ! Of if(chemthermod.ge.2)
313         
314         !Ion chemistry
315         if(chemthermod.eq.3) then
316            g_co2plus=igcm_co2plus
317            if (g_co2plus.eq.0 .and. chemthermod.ge.3) then
318               write(*,*) "chemthermos: Error; no CO2+ tracer !!!"
319               write(*,*) "CO2+ must be present if chemthermod=3"
320               stop
321            else
322               nesptherm=nesptherm+1 
323            endif
324            g_oplus=igcm_oplus
325            if (g_oplus.eq.0 .and. chemthermod.ge.3) then
326               write(*,*) "chemthermos: Error; no O+ tracer !!!"
327               write(*,*) "O+ must be present if chemthermod=3"
328               stop
329            else
330               nesptherm=nesptherm+1 
331            endif
332            g_o2plus=igcm_o2plus
333            if (g_o2plus.eq.0 .and. chemthermod.ge.3) then
334               write(*,*) "chemthermos: Error; no O2+ tracer !!!"
335               write(*,*) "O2+ must be present if chemthermod=3"
336               stop
337            else
338               nesptherm=nesptherm+1 
339            endif
340            g_coplus=igcm_coplus
341            if (g_coplus.eq.0) then
342               write(*,*) "chemthermos: Error; no CO+ tracer !!!"
343               write(*,*) "CO+ must be present if chemthermod=3"
344               stop
345            else
346               nesptherm=nesptherm+1 
347            endif
348            g_cplus=igcm_cplus
349            if (g_cplus.eq.0) then
350               write(*,*) "chemthermos: Error; no C+ tracer !!!"
351               write(*,*) "C+ must be present if chemthermod=3"
352               stop
353            else
354               nesptherm=nesptherm+1 
355            endif
356            g_nplus=igcm_nplus
357            if (g_nplus.eq.0) then
358               write(*,*) "chemthermos: Error; no N+ tracer !!!"
359               write(*,*) "N+ must be present if chemthermod=3"
360               stop
361            else
362               nesptherm=nesptherm+1 
363            endif
364            g_noplus=igcm_noplus
365            if (g_noplus.eq.0) then
366               write(*,*) "chemthermos: Error; no NO+ tracer !!!"
367               write(*,*) "NO+ must be present if chemthermod=3"
368               stop
369            else
370               nesptherm=nesptherm+1 
371            endif
372            g_n2plus=igcm_n2plus
373            if (g_n2plus.eq.0) then
374               write(*,*) "chemthermos: Error; no N2+ tracer !!!"
375               write(*,*) "N2+ must be present if chemthermod=3"
376               stop
377            else
378               nesptherm=nesptherm+1 
379            endif
380            g_hplus=igcm_hplus
381            if (g_hplus.eq.0) then
382               write(*,*) "chemthermos: Error; no H+ tracer !!!"
383               write(*,*) "H+ must be present if chemthermod=3"
384               stop
385            else
386               nesptherm=nesptherm+1 
387            endif
388            g_hco2plus=igcm_hco2plus
389            if (g_hco2plus.eq.0 .and. chemthermod.ge.3) then
390               write(*,*) "chemthermos: Error; no HCO2+ tracer !!!"
391               write(*,*) "HCO2+ must be present if chemthermod=3"
392               stop
393            else
394               nesptherm=nesptherm+1 
395            endif
396            g_elec=igcm_elec
397            if (g_elec.eq.0) then
398               write(*,*) "chemthermos: Error; no e- tracer !!!"
399               write(*,*) "e- must be present if chemthermod=3"
400               stop
401            else
402               nesptherm=nesptherm+1 
403            endif
404         endif
405
406         !Check if number of species is as expected
407         select case(chemthermod)
408         case(0)
409            if(nesptherm.ne.11) then
410               write(*,*)"chemthermos: Error:"
411               write(*,*)"Number of tracers not correct"
412               write(*,*)"There are",nesptherm," tracers"
413               write(*,*)"There should be 11 tracers"
414               stop
415            else
416               write(*,*)"chemthermos:",nesptherm," species"
417            endif
418         case(1)           
419            if(nesptherm.ne.12) then
420               write(*,*)"chemthermos: Error:"
421               write(*,*)"Number of tracers not correct"
422               write(*,*)"There are",nesptherm," tracers"
423               write(*,*)"There should be 12 tracers"
424               stop
425            else
426               write(*,*)"chemthermos:",nesptherm," species"
427            endif
428         case(2)
429            if(nesptherm.ne.17) then
430               write(*,*)"chemthermos: Error:"
431               write(*,*)"Number of tracers not correct"
432               write(*,*)"There are",nesptherm," tracers"
433               write(*,*)"There should be 17 tracers"
434               stop
435            else
436               write(*,*)"chemthermos:",nesptherm," species"
437            endif
438         case(3)
439            if(nesptherm.ne.28) then
440               write(*,*)"chemthermos: Error:"
441               write(*,*)"Number of tracers not correct"
442               write(*,*)"There are",nesptherm," tracers"
443               write(*,*)"There should be 28 tracers"
444               stop
445            else
446               write(*,*)"chemthermos:",nesptherm," species"
447            endif
448         end select
449         firstcall= .false.
450         write(*,*)'chemthermos: chemthermod=',chemthermod
451      endif                     ! of if (firstcall)
452     
453!cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
454
455      !Allocate density vector
[1266]456      allocate(rm(nlayer,nesptherm))
[635]457
[1266]458      do l=1,nlayer
[635]459        rm(l,i_co2)     = max(zycol(l,g_co2)*zdens(l),1.e-30)
460        rm(l,i_co)      = max(zycol(l,g_co)*zdens(l),1.e-30)
461        rm(l,i_o)       = max(zycol(l,g_o)*zdens(l),1.e-30)
462        rm(l,i_o1d)     = max(zycol(l,g_o1d)*zdens(l),1.e-30)
463        rm(l,i_o2)      = max(zycol(l,g_o2)*zdens(l),1.e-30)
464        rm(l,i_h)       = max(zycol(l,g_h)*zdens(l),1.e-30)
465        rm(l,i_h2)      = max(zycol(l,g_h2)*zdens(l),1.e-30)
466        rm(l,i_oh)      = max(zycol(l,g_oh)*zdens(l),1.e-30)
467        rm(l,i_ho2)     = max(zycol(l,g_ho2)*zdens(l),1.e-30)
468        rm(l,i_h2o2)    = max(zycol(l,g_h2o2)*zdens(l),1.e-30)
469        rm(l,i_h2o)     = max(zycol(l,g_h2o_vap)*zdens(l),1.e-30)
470        if(chemthermod.ge.1)  &
471             rm(l,i_o3)      = max(zycol(l,g_o3)*zdens(l),1.e-30)
472        if(chemthermod.ge.2) then
473           rm(l,i_n2)      = max(zycol(l,g_n2)*zdens(l),1.e-30)
474           rm(l,i_n)       = max(zycol(l,g_n)*zdens(l),1.e-30)
475           rm(l,i_no)      = max(zycol(l,g_no)*zdens(l),1.e-30)
476           rm(l,i_co)      = max(zycol(l,g_co)*zdens(l),1.e-30)
477           rm(l,i_n2d)     = max(zycol(l,g_n2d)*zdens(l),1.e-30)
478           rm(l,i_no2)     = max(zycol(l,g_no2)*zdens(l),1.e-30)
479        endif
480        if(chemthermod.eq.3) then
481           rm(l,i_co2plus) = max(zycol(l,g_co2plus)*zdens(l),1.e-30)
482           rm(l,i_oplus)   = max(zycol(l,g_oplus)*zdens(l),1.e-30)
483           rm(l,i_o2plus)  = max(zycol(l,g_o2plus)*zdens(l),1.e-30)
484           rm(l,i_coplus)  = max(zycol(l,g_coplus)*zdens(l),1.e-30)
485           rm(l,i_cplus)   = max(zycol(l,g_cplus)*zdens(l),1.e-30)
486           rm(l,i_nplus)   = max(zycol(l,g_nplus)*zdens(l),1.e-30)
487           rm(l,i_noplus)  = max(zycol(l,g_noplus)*zdens(l),1.e-30)
488           rm(l,i_n2plus)  = max(zycol(l,g_n2plus)*zdens(l),1.e-30)
489           rm(l,i_hplus)   = max(zycol(l,g_hplus)*zdens(l),1.e-30)
490           rm(l,i_hco2plus)= max(zycol(l,g_hco2plus)*zdens(l),1.e-30)
491           rm(l,i_elec)    = max(zycol(l,g_elec)*zdens(l),1.e-30)
492        endif
493        if(chemthermod.gt.3.or.chemthermod.lt.0) then
494           write(*,*)'chemthermos: bad value for chemthermod. Stop'
495           stop
496        endif
497      enddo
498
499      !Solar flux calculation
500
501      !Photoabsorption coefficients
[1684]502      call jthermcalc_e107(ig,nlayer,chemthermod,rm,nesptherm,ztemp,zlocal,zenit,zday)
[635]503
504      !Chemistry
505      call paramfoto_compact  &
[1266]506           (ig,nlayer,chemthermod,lswitch,ztemp,ptimestep,zenit,zlocal,rm,nesptherm)
[635]507
508      !Concentrations back to GCM
[1266]509      do l=lswitch,nlayer
[635]510        zycol(l,g_co2)     = max(rm(l,i_co2)     / zdens(l) , 1.e-30)
511        zycol(l,g_co)      = max(rm(l,i_co)      / zdens(l) , 1.e-30)
512        zycol(l,g_o2)      = max(rm(l,i_o2)      / zdens(l) , 1.e-30)
513        zycol(l,g_h2)      = max(rm(l,i_h2)      / zdens(l) , 1.e-30)
514        zycol(l,g_h)       = max(rm(l,i_h)       / zdens(l) , 1.e-30)
515        zycol(l,g_oh)      = max(rm(l,i_oh)      / zdens(l) , 1.e-30)
516        zycol(l,g_ho2)     = max(rm(l,i_ho2)     / zdens(l) , 1.e-30)
517        zycol(l,g_h2o_vap) = max(rm(l,i_h2o)     / zdens(l) , 1.e-30)
518        zycol(l,g_h2o2)    = max(rm(l,i_h2o2)    / zdens(l) , 1.e-30)
519        zycol(l,g_o1d)     = max(rm(l,i_o1d)     / zdens(l) , 1.e-30)
520        zycol(l,g_o)       = max(rm(l,i_o)       / zdens(l) , 1.e-30)
521        if(chemthermod.ge.1) &
522             zycol(l,g_o3) = max(rm(l,i_o3)      / zdens(l) , 1.e-30)
523        if(chemthermod.ge.2) then
524           zycol(l,g_n2)      = max(rm(l,i_n2)      / zdens(l) , 1.e-30)
525           zycol(l,g_n)       = max(rm(l,i_n)       / zdens(l) , 1.e-30)
526           zycol(l,g_no)      = max(rm(l,i_no)      / zdens(l) , 1.e-30)
527           zycol(l,g_no2)     = max(rm(l,i_no2)     / zdens(l) , 1.e-30)
528           zycol(l,g_n2d)     = max(rm(l,i_n2d)     / zdens(l) , 1.e-30)
529        endif
530        if(chemthermod.ge.3) then
531           zycol(l,g_co2plus) = max(rm(l,i_co2plus) / zdens(l) , 1.e-30)
532           zycol(l,g_oplus)   = max(rm(l,i_oplus)   / zdens(l) , 1.e-30)
533           zycol(l,g_o2plus)  = max(rm(l,i_o2plus)  / zdens(l) , 1.e-30)
534           zycol(l,g_coplus)  = max(rm(l,i_coplus)  / zdens(l) , 1.e-30)
535           zycol(l,g_cplus)   = max(rm(l,i_cplus)   / zdens(l) , 1.e-30)
536           zycol(l,g_nplus)   = max(rm(l,i_nplus)   / zdens(l) , 1.e-30)
537           zycol(l,g_noplus)  = max(rm(l,i_noplus)  / zdens(l) , 1.e-30)
538           zycol(l,g_n2plus)  = max(rm(l,i_n2plus)  / zdens(l) , 1.e-30)
539           zycol(l,g_hplus)   = max(rm(l,i_hplus)   / zdens(l) , 1.e-30)
540           zycol(l,g_hco2plus)= max(rm(l,i_hco2plus)/ zdens(l) , 1.e-30)
541           zycol(l,g_elec)    = max(rm(l,i_elec)    / zdens(l) , 1.e-30)
542        endif
[1888]543        em_no(l)=Pno(l,45)
544        em_o2(l)=Po2(l,10)*0.75
[635]545      enddo                     !nlayer
546
547      !Deallocations
548      deallocate(rm)
549
[3464]550      end subroutine chemthermos
[635]551
[3464]552  END MODULE chemthermos_mod
[635]553
554
555
Note: See TracBrowser for help on using the repository browser.