source: trunk/LMDZ.GENERIC/libf/aeronostd/concentrations.F @ 2537

Last change on this file since 2537 was 1796, checked in by bclmd, 7 years ago

Adding photochemistry to LMDZ Generic

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