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

Last change on this file since 3026 was 2615, checked in by romain.vande, 3 years ago

LMDZ_MARS RV : Open_MP;
Put all the "save" variables as "!$OMP THREADPRIVATE" in aeronomars

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