source: trunk/mars/libf/aeronomars/concentrations.F @ 86

Last change on this file since 86 was 38, checked in by emillour, 15 years ago

Ajout du modè Martien (mon LMDZ.MARS.BETA, du 28/01/2011) dans le rértoire mars, pour pouvoir suivre plus facilement les modifs.
EM
File size: 7.9 KB
RevLine 
[38]1      SUBROUTINE concentrations(pplay,pt,pdt,pq,pdq,ptimestep)
2                                             
3       IMPLICIT NONE
4c=======================================================================
5
6c CALCULATION OF MEAN MOLECULAR MASS, Cp, Akk and R
7c
8c mmean(ngridmx,nlayermx)                       amu
9c cpnew(ngridmx,nlayermx)                       J/kg/K
10c rnew(ngridmx,nlayermx)                        J/kg/K
11c akknew(ngridmx,nlayermx)                      coefficient of thermal concduction
12c
13c=======================================================================
14c    0.  Declarations :
15c    ------------------
16c
17#include "dimensions.h"
18#include "dimphys.h"
19#include "comcstfi.h"
20#include "callkeys.h"
21#include "comdiurn.h"
22#include "chimiedata.h"
23#include "tracer.h"
24#include "conc.h"
25
26c-----------------------------------------------------------------------
27c    Input/Output
28c    ------------
29
30      REAL pplay(ngridmx,nlayermx)
31      REAL pt(ngridmx,nlayermx)
32      REAL pdt(ngridmx,nlayermx)
33      real pq(ngridmx,nlayermx,nqmx)
34      REAL pdq(ngridmx,nlayermx,nqmx)
35      REAL ptimestep
36
37c    Local variables :
38c    -----------------
39      INTEGER,SAVE :: ngrid,nlayer,nq
40      INTEGER iq,l,ig,ll,n,k
41      integer,save :: gind(ncomptot)
42      real ni(ncomptot)
43      real nt, ntot
44      real q2(ngridmx,nlayermx,ncomptot)
45      real zt(ngridmx,nlayermx)
46      real q2tot(ngridmx,nlayermx)
47      real,save :: aki(ncomptot)
48      real,save :: cpi(ncomptot)
49
50      logical,save :: firstcall=.true.
51cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
52c     tracer numbering for the thermal conduction and
53c     specific heat coefficients
54cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
55c
56      integer,parameter :: i_co2  = 1
57      integer,parameter :: i_co   = 2
58      integer,parameter :: i_o    = 3
59      integer,parameter :: i_o1d  = 4
60      integer,parameter :: i_o2   = 5
61      integer,parameter :: i_o3   = 6
62      integer,parameter :: i_h    = 7
63      integer,parameter :: i_h2   = 8
64      integer,parameter :: i_oh   = 9
65      integer,parameter :: i_ho2  = 10
66      integer,parameter :: i_h2o2 = 11
67      integer,parameter :: i_n2   = 12
68      integer,parameter :: i_ar   = 13
69      integer,parameter :: i_h2o  = 14
70
71! Tracer indexes in the GCM:
72      integer,save :: g_co2=0
73      integer,save :: g_co=0
74      integer,save :: g_o=0
75      integer,save :: g_o1d=0
76      integer,save :: g_o2=0
77      integer,save :: g_o3=0
78      integer,save :: g_h=0
79      integer,save :: g_h2=0
80      integer,save :: g_oh=0
81      integer,save :: g_ho2=0
82      integer,save :: g_h2o2=0
83      integer,save :: g_n2=0
84      integer,save :: g_ar=0
85      integer,save :: g_h2o_vap=0
86
87! Initializations at first call
88      if (firstcall) then
89        ! identify the indexes of the tracers we'll need
90        g_co2=igcm_co2
91        if (g_co2.eq.0) then
92          write(*,*) "concentrations: Error; no CO2 tracer !!!"
93          stop
94        endif
95        g_co=igcm_co
96        if (g_co.eq.0) then
97          write(*,*) "concentrations: Error; no CO tracer !!!"
98          stop
99        endif
100        g_o=igcm_o
101        if (g_o.eq.0) then
102          write(*,*) "concentrations: Error; no O tracer !!!"
103          stop
104        endif
105        g_o1d=igcm_o1d
106        if (g_o1d.eq.0) then
107          write(*,*) "concentrations: Error; no O1D tracer !!!"
108          stop
109        endif
110        g_o2=igcm_o2
111        if (g_o2.eq.0) then
112          write(*,*) "concentrations: Error; no O2 tracer !!!"
113          stop
114        endif
115        g_o3=igcm_o3
116        if (g_o3.eq.0) then
117          write(*,*) "concentrations: Error; no O3 tracer !!!"
118          stop
119        endif
120        g_h=igcm_h
121        if (g_h.eq.0) then
122          write(*,*) "concentrations: Error; no H tracer !!!"
123          stop
124        endif
125        g_h2=igcm_h2
126        if (g_h2.eq.0) then
127          write(*,*) "concentrations: Error; no H2 tracer !!!"
128          stop
129        endif
130        g_oh=igcm_oh
131        if (g_oh.eq.0) then
132          write(*,*) "concentrations: Error; no OH tracer !!!"
133          stop
134        endif
135        g_ho2=igcm_ho2
136        if (g_ho2.eq.0) then
137          write(*,*) "concentrations: Error; no HO2 tracer !!!"
138          stop
139        endif
140        g_h2o2=igcm_h2o2
141        if (g_h2o2.eq.0) then
142          write(*,*) "concentrations: Error; no H2O2 tracer !!!"
143          stop
144        endif
145        g_n2=igcm_n2
146        if (g_n2.eq.0) then
147          write(*,*) "concentrations: Error; no N2 tracer !!!"
148          stop
149        endif
150        g_ar=igcm_ar
151        if (g_ar.eq.0) then
152          write(*,*) "concentrations: Error; no AR tracer !!!"
153          stop
154        endif
155        g_h2o_vap=igcm_h2o_vap
156        if (g_h2o_vap.eq.0) then
157          write(*,*) "concentrations: Error; no water vapor tracer !!!"
158          stop
159        endif
160
161cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
162c     fill local array of tracer indexes
163cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
164c
165        gind(i_co2)   =  g_co2                  ! co2
166        gind(i_co)    =  g_co                   ! co
167        gind(i_o)     =  g_o                    ! o
168        gind(i_o1d)   =  g_o1d                  ! o1d
169        gind(i_o2)    =  g_o2                   ! o2
170        gind(i_o3)    =  g_o3                   ! o3
171        gind(i_h)     =  g_h                    ! h
172        gind(i_h2)    =  g_h2                   ! h2
173        gind(i_oh)    =  g_oh                   ! oh
174        gind(i_ho2)   =  g_ho2                  ! ho2
175        gind(i_h2o2)  =  g_h2o2                 ! h2o2
176        gind(i_n2)    =  g_n2                   ! n2
177        gind(i_ar)    =  g_ar                   ! ar
178        gind(i_h2o)   =  g_h2o_vap              ! h2o
179
180cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
181c     Thermal conductivity and specific heat coefficients
182cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
183c
184        aki(i_co2)  = 3.072e-4
185        aki(i_co)   = 4.87e-4
186        aki(i_o)    = 7.59e-4
187        aki(i_o1d)  = 7.59e-4           !?
188        aki(i_o2)   = 5.68e-4
189        aki(i_o3)   = 3.00e-4           !?
190        aki(i_h)    = 0.0
191        aki(i_h2)   = 36.314e-4
192        aki(i_oh)   = 7.00e-4           !?
193        aki(i_ho2)  = 0.0
194        aki(i_h2o2) = 0.0
195        aki(i_n2)   = 5.6e-4
196        aki(i_ar)   = 0.0                 !?
197        aki(i_h2o)  = 0.0
198
199        cpi(i_co2)   = 0.834e3
200        cpi(i_co)    = 1.034e3
201        cpi(i_o)     = 1.3e3
202        cpi(i_o1d)   = 1.3e3            !?
203        cpi(i_o2)    = 0.9194e3
204        cpi(i_o3)    = 0.800e3  !?
205        cpi(i_h)     = 20.780e3
206        cpi(i_h2)    = 14.266e3
207        cpi(i_oh)    = 1.045e3
208        cpi(i_ho2)   = 1.065e3  !?
209        cpi(i_h2o2)  = 1.000e3  !?
210        cpi(i_n2)    = 1.034e3
211        cpi(i_ar)    = 1.000e3    !?
212        cpi(i_h2o)   = 1.870e3
213c
214cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
215
216        nlayer=nlayermx
217        ngrid=ngridmx
218        nq=nqmx
219
220        firstcall=.false.
221      endif ! of if (firstcall)
222
223      DO l=1,nlayer
224        DO ig=1,ngrid
225          DO n=1,ncomptot
226            q2(ig,l,n)=max(1.e-30,
227     .         pq(ig,l,gind(n))+pdq(ig,l,gind(n))*ptimestep)
228          ENDDO
229          zt(ig,l)=pt(ig,l) +pdt(ig,l)*ptimestep
230        ENDDO
231      ENDDO
232 
233      do l=1,nlayermx
234        do ig=1,ngridmx
235          ntot=pplay(ig,l)/(1.381e-23*zt(ig,l))*1.e-6  ! in #/cm3
236          cpnew(ig,l)=0.
237          akknew(ig,l)=0.
238          mmean(ig,l)=0.
239          q2tot(ig,l)=0.
240          nt=0.
241          do n=1,ncomptot
242            ni(n)=0.0
243            do k=1,ncomptot
244              if(k.ne.n) ni(n)=ni(n)+q2(ig,l,k)/mmol(gind(k))
245            enddo
246            ni(n)=ntot/(1.+mmol(gind(n))/q2(ig,l,n)*ni(n))
247            cpnew(ig,l)=cpnew(ig,l)+ni(n)*cpi(n)
248            akknew(ig,l)=akknew(ig,l)+ni(n)*aki(n)
249            mmean(ig,l)=mmean(ig,l)+q2(ig,l,n)/mmol(gind(n))
250            q2tot(ig,l)=q2tot(ig,l)+q2(ig,l,n)
251c        if(ig.eq.1.and.l.eq.1) write(*,*)'q2tot(ig,l)',n,q2tot(ig,l)
252            if(cpi(n) .ne. 0.0) nt=nt+ni(n)
253          enddo
254
255c     print*,"concentrations rep 3",l,ig,nt,mmean(ig,l),zt(ig+1,l)
256
257          cpnew(ig,l)=cpnew(ig,l)/nt
258          akknew(ig,l)=akknew(ig,l)/nt
259          mmean(ig,l)=1/mmean(ig,l)    ! in amu
260          rnew(ig,l)=8.314/mmean(ig,l)*1.e3     ! J/kg/K               
261        enddo
262c        print*,l,mmean(1,l),cpnew(1,l),rnew(1,l)
263c         write(228,*),l,pplay(1,l),ntot
264      enddo
265      return
266      end
Note: See TracBrowser for help on using the repository browser.