source: trunk/LMDZ.MARS/libf/aeronomars/concentrations.F @ 2599

Last change on this file since 2599 was 2158, checked in by emillour, 5 years ago

Mars GCM:

  • Updated chemical core to include ionospheric chemistry

FGG

File size: 9.3 KB
Line 
1      SUBROUTINE concentrations(ngrid,nlayer,nq,
2     &                          pplay,pt,pdt,pq,pdq,ptimestep)
3                                             
4      use tracer_mod, only: igcm_co2, igcm_co, igcm_o, igcm_o1d,
5     &                      igcm_o2, igcm_o3, igcm_h, igcm_h2,
6     &                      igcm_oh, igcm_ho2, igcm_n2, igcm_ar,
7     &                      igcm_h2o_vap, igcm_n, igcm_no, igcm_no2,
8     &                      igcm_n2d, igcm_co2plus, igcm_oplus,
9     &                      igcm_o2plus, igcm_coplus, igcm_cplus,
10     &                      igcm_nplus, igcm_noplus, igcm_n2plus,
11     &                      igcm_hplus, igcm_hco2plus, mmol,
12     &                      igcm_he, igcm_elec
13      use conc_mod, only: mmean, Akknew, rnew, cpnew
14      implicit none
15
16!=======================================================================
17! CALCULATION OF MEAN MOLECULAR MASS, Cp, Akk and R
18!
19! mmean(ngrid,nlayer)   amu
20! cpnew(ngrid,nlayer)   J/kg/K
21! rnew(ngrid,nlayer)    J/kg/K
22! akknew(ngrid,nlayer)  coefficient of thermal concduction
23!
24! version: April 2012 - Franck Lefevre
25!=======================================================================
26
27!     declarations
28 
29#include "callkeys.h"
30
31!     input/output
32
33      integer,intent(in) :: ngrid ! number of atmospheric columns
34      integer,intent(in) :: nlayer ! number of atmospheric layers
35      integer,intent(in) :: nq ! number of tracers
36      real,intent(in) :: pplay(ngrid,nlayer)
37      real,intent(in) :: pt(ngrid,nlayer)
38      real,intent(in) :: pdt(ngrid,nlayer)
39      real,intent(in) :: pq(ngrid,nlayer,nq)
40      real,intent(in) :: pdq(ngrid,nlayer,nq)
41      real,intent(in) :: ptimestep
42
43!     local variables
44
45      integer       :: i, l, ig, iq
46      integer, save :: nbq
47      integer,allocatable,save :: niq(:)
48      real          :: ni(nq), ntot
49      real          :: zq(ngrid, nlayer, nq)
50      real          :: zt(ngrid, nlayer)
51      real,allocatable,save    :: aki(:)
52      real,allocatable,save    :: cpi(:)
53
54      logical, save :: firstcall = .true.
55
56      if (firstcall) then
57
58         ! allocate local saved arrays:
59         allocate(aki(nq))
60         allocate(cpi(nq))
61         allocate(niq(nq))
62!        find index of chemical tracers to use
63!        initialize thermal conductivity and specific heat coefficients
64!        !? values are estimated
65
66         nbq = 0 ! to count number of tracers used in this subroutine
67
68         if (igcm_co2 /= 0) then
69            nbq = nbq + 1
70            niq(nbq) = igcm_co2
71            aki(nbq) = 3.072e-4
72            cpi(nbq) = 0.834e3
73         end if
74         if (igcm_co /= 0) then
75            nbq = nbq + 1
76            niq(nbq) = igcm_co
77            aki(nbq) = 4.87e-4
78            cpi(nbq) = 1.034e3
79         end if
80         if (igcm_o /= 0) then
81            nbq = nbq + 1
82            niq(nbq) = igcm_o
83            aki(nbq) = 7.59e-4
84            cpi(nbq) = 1.3e3
85         end if
86         if (igcm_o1d /= 0) then
87            nbq = nbq + 1
88            niq(nbq) = igcm_o1d
89            aki(nbq) = 7.59e-4  !?
90            cpi(nbq) = 1.3e3    !?
91         end if
92         if (igcm_o2 /= 0) then
93            nbq = nbq + 1
94            niq(nbq) = igcm_o2
95            aki(nbq) = 5.68e-4
96            cpi(nbq) = 0.9194e3
97         end if
98         if (igcm_o3 /= 0) then
99            nbq = nbq + 1
100            niq(nbq) = igcm_o3
101            aki(nbq) = 3.00e-4  !?
102            cpi(nbq) = 0.800e3  !?
103         end if
104         if (igcm_h /= 0) then
105            nbq = nbq + 1
106            niq(nbq) = igcm_h
107            aki(nbq) = 0.0
108            cpi(nbq) = 20.780e3
109         end if
110         if (igcm_h2 /= 0) then
111            nbq = nbq + 1
112            niq(nbq) = igcm_h2
113            aki(nbq) = 36.314e-4
114            cpi(nbq) = 14.266e3
115         end if
116         if (igcm_oh /= 0) then
117            nbq = nbq + 1
118            niq(nbq) = igcm_oh
119            aki(nbq)  = 7.00e-4 !?
120            cpi(nbq)  = 1.045e3
121         end if
122         if (igcm_ho2 /= 0) then
123            nbq = nbq + 1
124            niq(nbq) = igcm_ho2
125            aki(nbq) = 0.0
126            cpi(nbq) = 1.065e3  !?
127         end if
128         if (igcm_n2 /= 0) then
129            nbq = nbq + 1
130            niq(nbq) = igcm_n2
131            aki(nbq) = 5.6e-4
132            cpi(nbq) = 1.034e3
133         end if
134         if (igcm_ar /= 0) then
135            nbq = nbq + 1
136            niq(nbq) = igcm_ar
137            aki(nbq) = 3.4e-4
138            cpi(nbq) = 5.2e2
139            ! aki/cpi from Leslie A. Guildner,
140            ! JOURNAL OF RESEARCH of the National Bureau of Standards-
141            ! A. Physics and Chemistry
142            ! Vol. 79A, No.2, March-April 1975
143         end if
144         if (igcm_h2o_vap /= 0) then
145            nbq = nbq + 1
146            niq(nbq) = igcm_h2o_vap
147            aki(nbq) = 0.0
148            cpi(nbq) = 1.870e3
149         end if
150         if (igcm_n /= 0) then
151            nbq = nbq + 1
152            niq(nbq) = igcm_n
153            aki(nbq) = 0.0
154            cpi(nbq) = 0.0
155         endif
156         if(igcm_no /= 0) then
157            nbq = nbq + 1
158            niq(nbq) = igcm_no
159            aki(nbq) = 0.0
160            cpi(nbq) = 0.0
161         endif
162         if(igcm_no2 /= 0) then
163            nbq = nbq + 1
164            niq(nbq) = igcm_no2
165            aki(nbq) = 0.0
166            cpi(nbq) = 0.0
167         endif
168         if(igcm_n2d /= 0) then
169            nbq = nbq + 1
170            niq(nbq) = igcm_n2d
171            aki(nbq) = 0.0
172            cpi(nbq) = 0.0
173         endif
174         if(igcm_he /= 0) then
175            nbq = nbq + 1
176            niq(nbq) = igcm_he
177            aki(nbq) = 30.e-4
178            cpi(nbq) = 5.2e3
179            ! aki/cpi from Leslie A. Guildner,
180            ! JOURNAL OF RESEARCH of the National Bureau of Standards-
181            ! A. Physics and Chemistry
182            ! Vol. 79A, No.2, March-April 1975
183         endif
184         if(igcm_co2plus /= 0) then
185            nbq = nbq + 1
186            niq(nbq) = igcm_co2plus
187            aki(nbq) = 0.0
188            cpi(nbq) = 0.0
189         endif
190         if(igcm_oplus /= 0) then
191            nbq = nbq + 1
192            niq(nbq) = igcm_oplus
193            aki(nbq) = 0.0
194            cpi(nbq) = 0.0
195         endif
196         if(igcm_o2plus /= 0) then
197            nbq = nbq + 1
198            niq(nbq) = igcm_o2plus
199            aki(nbq) = 0.0
200            cpi(nbq) = 0.0
201         endif
202         if(igcm_coplus /= 0) then
203            nbq = nbq + 1
204            niq(nbq) = igcm_coplus
205            aki(nbq) = 0.0
206            cpi(nbq) = 0.0
207         endif
208         if(igcm_cplus /= 0) then
209            nbq = nbq + 1
210            niq(nbq) = igcm_cplus
211            aki(nbq) = 0.0
212            cpi(nbq) = 0.0
213         endif
214         if(igcm_nplus /= 0) then
215            nbq = nbq + 1
216            niq(nbq) = igcm_nplus
217            aki(nbq) = 0.0
218            cpi(nbq) = 0.0
219         endif
220         if(igcm_noplus /= 0) then
221            nbq = nbq + 1
222            niq(nbq) = igcm_noplus
223            aki(nbq) = 0.0
224            cpi(nbq) = 0.0
225         endif
226         if(igcm_n2plus /= 0) then
227            nbq = nbq + 1
228            niq(nbq) = igcm_n2plus
229            aki(nbq) = 0.0
230            cpi(nbq) = 0.0
231         endif
232         if(igcm_hplus /= 0) then
233            nbq = nbq + 1
234            niq(nbq) = igcm_hplus
235            aki(nbq) = 0.0
236            cpi(nbq) = 0.0
237         endif
238         if(igcm_hco2plus /= 0) then
239            nbq = nbq + 1
240            niq(nbq) = igcm_hco2plus
241            aki(nbq) = 0.0
242            cpi(nbq) = 0.0
243         endif
244         if(igcm_elec /= 0) then
245            nbq = nbq + 1
246            niq(nbq) = igcm_elec
247            aki(nbq) = 0.0
248            cpi(nbq) = 0.0
249         endif
250         
251         ! tell the world about it:
252         write(*,*) "concentrations: firstcall, nbq=",nbq
253         write(*,*) "  niq(1:nbq)=",niq(1:nbq)
254         write(*,*) "  aki(1:nbq)=",aki(1:nbq)
255         write(*,*) "  cpi(1:nbq)=",cpi(1:nbq)
256
257         firstcall = .false.
258
259      end if ! if (firstcall)
260
261!     update temperature
262
263      do l = 1,nlayer
264         do ig = 1,ngrid
265            zt(ig,l) = pt(ig,l) + pdt(ig,l)*ptimestep
266         end do
267      end do
268
269!     update tracers
270
271      do l = 1,nlayer
272         do ig = 1,ngrid
273            do i = 1,nbq
274               iq = niq(i)
275               zq(ig,l,iq) = max(1.e-30, pq(ig,l,iq)
276     $                                 + pdq(ig,l,iq)*ptimestep)
277            end do
278         end do
279      end do
280
281!     mmean : mean molecular mass
282!     rnew  : specific gas constant
283
284      mmean(:,:)  = 0.
285
286      do l = 1,nlayer
287         do ig = 1,ngrid
288            do i = 1,nbq
289               iq = niq(i)
290               mmean(ig,l) = mmean(ig,l) + zq(ig,l,iq)/mmol(iq)
291            end do
292            mmean(ig,l) = 1./mmean(ig,l)
293            rnew(ig,l) = 8.314/mmean(ig,l)*1.e3     ! J/kg/K           
294         end do
295      end do
296
297!     cpnew  : specicic heat
298!     akknew : thermal conductivity cofficient
299     
300      cpnew(:,:)  = 0.
301      akknew(:,:) = 0.
302
303      do l = 1,nlayer
304         do ig = 1,ngrid
305            ntot = pplay(ig,l)/(1.381e-23*zt(ig,l))*1.e-6  ! in #/cm3
306            do i = 1,nbq
307               iq = niq(i)
308               ni(iq) = ntot*zq(ig,l,iq)*mmean(ig,l)/mmol(iq)
309               cpnew(ig,l) = cpnew(ig,l) + ni(iq)*cpi(i)
310               akknew(ig,l) = akknew(ig,l) + ni(iq)*aki(i)
311            end do
312            cpnew(ig,l) = cpnew(ig,l)/ntot
313            akknew(ig,l) = akknew(ig,l)/ntot
314         end do
315!        print*, l, mmean(1,l), cpnew(1,l), rnew(1,l)
316      end do
317
318      return
319      end
Note: See TracBrowser for help on using the repository browser.