source: trunk/LMDZ.MARS/libf/phymars/nltecool.F @ 3759

Last change on this file since 3759 was 3726, checked in by emillour, 10 months ago

Mars PCM:
Turn "callkeys.h" into module "callkeys_mod.F90"
EM
File size: 18.9 KB
Line 
1      MODULE nltecool_mod
2     
3      IMPLICIT NONE
4
5      ! Escape functions and VMRs from tabulated values.
6      ! Origin: nlte_escape.dat file form Miguel L.
7      ! (Y. Wanherdrick, 09/2000)
8      integer,parameter :: np=68 ! # data points in tabulation
9     
10      real, save :: pnb(np) = [   ! Pressure (exponent) in tabulation
11     &    12.0000,    11.0000,    10.8000,
12     &    10.6000,   10.40000,   10.20000,
13     &   10.00000,    9.80000,    9.60000,
14     &    9.40000,    9.20000,    9.00000,
15     &    8.80000,    8.60000,    8.40000,
16     &    8.20000,    8.00000,    7.80000,
17     &    7.60000,    7.40000,    7.20000,
18     &    7.00000,    6.80000,    6.60000,
19     &    6.40000,    6.20000,    6.00000,
20     &    5.80000,    5.60000,    5.40000,
21     &    5.20000,    5.00000,    4.80000,
22     &    4.60000,    4.40000,    4.20000,
23     &    4.00000,    3.80000,    3.60000,
24     &    3.40000,    3.20000,    3.00000,
25     &    2.80000,    2.60000,    2.40000,
26     &    2.20000,    2.00000,    1.80000,
27     &    1.60000,    1.40000,    1.20000,
28     &    1.00000,   0.800000,   0.599999,
29     &   0.400000,   0.200000,  0.,
30     &  -0.200000,  -0.400001,  -0.600000,
31     &  -0.800000,   -1.00000,   -1.20000,
32     &   -1.40000,   -1.60000,   -1.80000,
33     &   -2.00000,   -3.00000 ] ! values modified in nltecool
34!$OMP THREADPRIVATE(pnb)
35
36      real,parameter :: ef1(np) = [   ! Esc.funct.#1, tabulated
37     &    4.58112E-04,    4.58112E-04,    4.58112E-04,
38     &    4.58112E-04,    4.58112E-04,    4.58112E-04,
39     &    4.58112E-04,    4.58112E-04,    4.58112E-04,
40     &    4.58112E-04,    4.58112E-04,    4.58112E-04,
41     &    4.58112E-04,    4.58112E-04,    4.58112E-04,
42     &    4.58112E-04,    4.58112E-04,    4.58112E-04,
43     &    4.58112E-04,    4.58112E-04,    4.58112E-04,
44     &    4.58112E-04,    4.61707E-04,    4.76886E-04,
45     &    4.95638E-04,    5.20935E-04,    5.55511E-04,
46     &    6.01219E-04,    6.63734E-04,    7.50691E-04,
47     &    8.63474E-04,    1.00900E-03,    1.19642E-03,
48     &    1.42690E-03,    1.71398E-03,    2.06663E-03,
49     &    2.48974E-03,    3.01578E-03,    3.64350E-03,
50     &    4.40323E-03,    5.32066E-03,    6.40456E-03,
51     &    7.72069E-03,    9.25684E-03,    1.10905E-02,
52     &    1.32374E-02,    1.57643E-02,    1.87388E-02,
53     &    2.22072E-02,    2.63099E-02,    3.10614E-02,
54     &    3.66948E-02,    4.32373E-02,    5.15022E-02,
55     &    6.21455E-02,    7.77212E-02,    9.92027E-02,
56     &   0.131155,   0.179470,   0.258913,
57     &   0.380549,   0.530450,   0.643180,
58     &   0.741061,   0.826336,   0.922787,
59     &   0.997203,    1.00000 ]
60     
61      real,parameter :: ef2(np) = [   ! Esc.funct.#2, tabulated
62     &    1.98992E-03,    1.98992E-03,    1.98992E-03,
63     &    1.98992E-03,    1.98992E-03,    1.98992E-03,
64     &    1.98992E-03,    1.98992E-03,    1.98992E-03,
65     &    1.98992E-03,    1.98992E-03,    2.01376E-03,
66     &    2.09450E-03,    2.22993E-03,    2.42056E-03,
67     &    2.68018E-03,    3.04398E-03,    3.43896E-03,
68     &    3.80282E-03,    4.20622E-03,    4.76121E-03,
69     &    8.01698E-03,    1.19947E-02,    1.69149E-02,
70     &    2.24497E-02,    2.85244E-02,    3.54813E-02,
71     &    4.39264E-02,    5.46248E-02,    6.75367E-02,
72     &    8.29931E-02,    1.01717E-01,   0.123422,
73     &   0.148468,   0.177096,   0.208816,
74     &   0.244003,   0.282013,   0.322559,
75     &   0.365542,   0.410518,   0.457384,
76     &   0.505358,   0.553627,   0.600472,
77     &   0.644807,   0.687185,   0.727429,
78     &   0.764734,   0.798562,   0.828699,
79     &   0.854797,   0.877717,   0.897874,
80     &   0.915258,   0.929904,   0.942381,
81     &   0.952906,   0.962173,   0.970191,
82     &   0.976437,   0.981501,   0.985406,
83     &   0.988560,   0.991111,   0.993653,
84     &   0.995561,    1.00000 ]
85     
86      real,parameter :: co2vmr(np) = [   ! CO2 VMR tabulated
87     &   0.950000,   0.950000,   0.950000,
88     &   0.950000,   0.950000,   0.950000,
89     &   0.950000,   0.950000,   0.950000,
90     &   0.950000,   0.950000,   0.950000,
91     &   0.950000,   0.950000,   0.950000,
92     &   0.950000,   0.950000,   0.950000,
93     &   0.950000,   0.950000,   0.950000,
94     &   0.950000,   0.950000,   0.950000,
95     &   0.950000,   0.950000,   0.950000,
96     &   0.950000,   0.950000,   0.950000,
97     &   0.950000,   0.950000,   0.950000,
98     &   0.950000,   0.950000,   0.950000,
99     &   0.950000,   0.950000,   0.950000,
100     &   0.950000,   0.950000,   0.950000,
101     &   0.950000,   0.950000,   0.950000,
102     &   0.950000,   0.950000,   0.950000,
103     &   0.950000,   0.950000,   0.950000,
104     &   0.950000,   0.950000,   0.950000,
105     &   0.950000,   0.950000,   0.950000,
106     &   0.950000,   0.950000,   0.950000,
107     &   0.949619,   0.947694,   0.945830,
108     &   0.944016,   0.940557,   0.937068,
109     &   0.932366,   0.893661 ]     
110
111      real,parameter :: o3pvmr(np) = [   ! O3p VMR tabulated
112     &    5.06756E-08,    9.16539E-07,    1.68217E-06,
113     &    3.00843E-06,    5.03151E-06,    8.07489E-06,
114     &    1.23137E-05,    1.79029E-05,    2.45308E-05,
115     &    3.27431E-05,    4.26692E-05,    5.44396E-05,
116     &    6.78865E-05,    8.33147E-05,    1.00148E-04,
117     &    1.18846E-04,    1.39681E-04,    1.64909E-04,
118     &    1.93617E-04,    2.25161E-04,    2.60834E-04,
119     &    3.01501E-04,    3.44953E-04,    3.91011E-04,
120     &    4.40377E-04,    4.90820E-04,    5.43200E-04,
121     &    5.95335E-04,    6.45420E-04,    6.93166E-04,
122     &    7.43729E-04,    7.93710E-04,    8.44394E-04,
123     &    8.94318E-04,    9.44732E-04,    9.94964E-04,
124     &    1.04901E-03,    1.10008E-03,    1.16302E-03,
125     &    1.22989E-03,    1.30026E-03,    1.37131E-03,
126     &    1.45556E-03,    1.55186E-03,    1.66328E-03,
127     &    1.77802E-03,    1.91546E-03,    2.07503E-03,
128     &    2.24903E-03,    2.47117E-03,    2.71728E-03,
129     &    2.99739E-03,    3.33582E-03,    3.73507E-03,
130     &    4.20819E-03,    4.76887E-03,    5.42558E-03,
131     &    6.20815E-03,    7.14473E-03,    8.28545E-03,
132     &    9.51779E-03,    1.08140E-02,    1.22359E-02,
133     &    1.36870E-02,    1.51495E-02,    1.67196E-02,
134     &    1.85485E-02,    3.36252E-02 ]
135
136      real, parameter :: n2covmr(np) = [   ! N2+CO VMR tabulated
137     &    2.71412E-02,    2.71464E-02,    2.71490E-02,
138     &    2.71523E-02,    2.71558E-02,    2.71617E-02,
139     &    2.71672E-02,    2.71749E-02,    2.71837E-02,
140     &    2.71943E-02,    2.72058E-02,    2.72189E-02,
141     &    2.72326E-02,    2.72483E-02,    2.72661E-02,
142     &    2.72848E-02,    2.73054E-02,    2.73279E-02,
143     &    2.73514E-02,    2.73775E-02,    2.74048E-02,
144     &    2.74345E-02,    2.74672E-02,    2.75021E-02,
145     &    2.75404E-02,    2.75826E-02,    2.76340E-02,
146     &    2.77013E-02,    2.78220E-02,    2.79707E-02,
147     &    2.81759E-02,    2.84339E-02,    2.87587E-02,
148     &    2.91600E-02,    2.96561E-02,    3.02558E-02,
149     &    3.09922E-02,    3.18062E-02,    3.27010E-02,
150     &    3.35635E-02,    3.44388E-02,    3.53327E-02,
151     &    3.62143E-02,    3.70941E-02,    3.79315E-02,
152     &    3.87626E-02,    3.95524E-02,    4.03071E-02,
153     &    4.10071E-02,    4.16229E-02,    4.21231E-02,
154     &    4.25167E-02,    4.27964E-02,    4.29773E-02,
155     &    4.30488E-02,    4.29638E-02,    4.28049E-02,
156     &    4.26788E-02,    4.26822E-02,    4.29426E-02,
157     &    4.34634E-02,    4.42559E-02,    4.53038E-02,
158     &    4.65879E-02,    4.80262E-02,    4.96303E-02,
159     &    5.14885E-02,    6.91651E-02 ]
160
161      CONTAINS
162c**************************************************************************
163c
164      subroutine nltecool(ngrid,nlayer,nq,pplay,pt,pq,dtnlte)
165c
166c This code was designed as a delivery for the "Martian Environment Models"
167c project ( ESA contract 11369/95/nl/jg CCN2 )
168c Computes non-LTE heating rates from CO2 emission at 15 um
169c in the Martian upper atmosphere.
170c Uses a simplified model consisting of two excited levels with two
171c emission bands, one of them stronger than the other, which correspond
172c to the behaviours of the 626 fundamental band and the isotopic fund.bands.
173c It uses a cool-to-space approximation with tabulated escape functions.
174c These escape functions have been precomputed for the strong and weak bands,
175c and are given as a function of pressure in separate files.
176c The output values are the heating rates (actually, cooling, since they
177c are always negative) for the two bands, i.e., the total cooling is the
178c sum of them.
179c Miguel A. Lopez Valverde
180c Instituto de Astrofisica de Andalucia (CSIC), Granada, Spain
181c
182c Version 1b.  See description above.  22-March-2000.
183c Adapted as a subroutine for use in GCM -- PLR/SRL 6/2000
184c Version 1c.  Inclusion of VMR in the tabulation of escape functions.
185c              Table contains now only 1 input file -- Miguel 11/Jul/2000
186c Version 1d  data contained in original input file "nlte_escape.dat"
187c now stored in include file "nltedata.h" Y.Wanherdrick 09/2000
188
189c       jul 2011 fgg   Modified to allow variable O
190c     
191c***************************************************************************
192
193      use tracer_mod, only: igcm_co2, igcm_co, igcm_o, igcm_n2, mmol
194      use conc_mod, only: mmean
195      use callkeys_mod, only: nltemodel
196      implicit none
197
198c Input and output variables
199c
200      integer,intent(in) :: ngrid               ! no. of horiz. gridpoints
201      integer,intent(in) :: nlayer              ! no. of atmospheric layers
202      integer,intent(in) :: nq                  ! no. of tracers
203      real,intent(in) :: pplay(ngrid,nlayer)    ! input pressure grid (Pa)
204      real,intent(in) :: pt(ngrid,nlayer)       ! input temperatures (K)
205      real,intent(in) :: pq(ngrid,nlayer,nq)    ! input mmrs (kg/kg_air)
206      real,intent(out) :: dtnlte(ngrid,nlayer)  ! output temp. tendencies (K/s)
207
208c
209c Standard atmosphere variables
210c
211      real        nt                              ! number density [cm-3]
212      real        co2(nlayer)                     !  "   of CO2
213      real        o3p(nlayer)                     !  "   of atomic oxygen
214      real        n2co(nlayer)                    !  "  of N2 + CO
215      real        pyy(nlayer)                     ! auxiliary pressure grid
216
217c
218c Vectors and indexes for the tabulation of escape functions and VMR
219c
220c                 np                              ! # data points in tabulation
221c                 pnb(np)                         ! Pressure in tabulation
222c                 ef1(np)                         ! Esc.funct.#1, tabulated
223c                 ef2(np)                         ! Esc.funct.#2, tabulated
224c                 co2vmr(np)                      ! CO2 VMR tabulated
225c                 o3pvmr(np)                      ! CO2 VMR tabulated
226c                 n2covmr(np)                     ! N2+CO VMR tabulated
227      real        escf1(nlayer)                   ! Esc.funct.#1, interpolated
228      real        escf2(nlayer)                   ! Esc.funct.#2, interpolated
229
230
231c
232c Local Constants
233c
234      real       nu1, nu2                         ! freq. of energy levels
235      real       imr1, imr2                       ! isotopic abundances
236      real       hplanck, gamma, vlight           ! physical constants
237      real       ee
238      real       rfvt                             ! collisional rate
239      real       rfvto3p                          !     "
240      real       rfvv                             !     "
241
242c
243c Local variables for the main loop
244c
245      real       n1, n2, co2t                     ! ground populations
246      real       l1, p1, p12                      ! prod & losses
247      real       l2, p2, p21
248      real       tt                               ! dummy variable
249      real       c1, c2                           ! molecular constants
250      real       ae1, ae2                         ! einstein spontaneous emission
251      real       a1, a2, a12, a21
252      real       pl1, pl2
253      real       el1, el2
254      real       hr1, hr2                         ! heating rate due to each band
255      real       hr(nlayer)                       ! total heating rate
256
257c
258c Indexes
259c
260      integer    i
261      integer    j,ii
262
263c
264c Rate coefficients
265c
266      real       k19xca, k19xcb
267      real       k19cap1, k19cap2
268      real       k19cbp1, k19cbp2
269      real       d19c, d19cp1, d19cp2
270      real       k20xc, k20cp1, k20cp2
271      real       k21xc, k21cp2
272
273      logical,save :: firstcall=.true.
274!$OMP THREADPRIVATE(firstcall)
275
276c
277c Data
278c
279      data       nu1, nu2, hplanck, gamma, vlight, ee/
280     1     667.38, 662.3734, 6.6261e-27, 1.191e-5, 3.e10, 1.438769/
281
282c*************************************************************************
283c       PROGRAM  STARTS
284c*************************************************************************
285
286      imr1 = 0.987
287      imr2 = 0.00408 + 0.0112
288      rfvt = 0.1
289      rfvto3p = 1.0
290      rfvv = 0.1
291
292      if(firstcall) then
293
294         do i=1,np
295            pnb(i)=1.0e-4*exp(pnb(i)) ! p into Pa
296         end do
297
298         firstcall = .false.
299
300      endif
301c
302c MAIN LOOP, for each gridpoint and altitude:
303c
304      do j=1,ngrid  ! loop over grid points
305c
306c set up local pressure grid
307c
308         do ii=1,nlayer
309            pyy(ii)=pplay(j,ii)
310         enddo
311!
312! Interpolate escape functions and VMR to the desired grid
313!
314         call interp1(escf2,pyy,nlayer,ef2,pnb,np)
315         call interp1(escf1,pyy,nlayer,ef1,pnb,np)
316         if(nltemodel.eq.0) then
317            call interp3(co2,o3p,n2co,pyy,nlayer,
318     &           co2vmr,o3pvmr,n2covmr,pnb,np)
319         endif
320       
321         do i=1,nlayer  ! loop over layers
322C
323C test if p lies outside range (p > 3.5 Pa)
324C changed to 1 Pa since transition will always be higher than this
325C
326            if(pyy(i) .gt. 1.0 .or. pyy(i) .lt. 4.0e-6) then
327               hr(i)=0.0
328               dtnlte(j,i)=0.0
329            else
330c
331c           if(pt(j,i).lt.1.0)print*,pt(j,i)
332               nt = pyy(i)/(1.381e-17*pt(j,i)) ! nt in cm-3
333               !Dynamic composition
334               if(nltemodel.eq.1) then
335                  co2(i)=pq(j,i,igcm_co2)*mmean(j,i)/mmol(igcm_co2)
336                  o3p(i)=pq(j,i,igcm_o)*mmean(j,i)/mmol(igcm_o)
337                  n2co(i)=pq(j,i,igcm_co)*mmean(j,i)/mmol(igcm_co) +
338     $                 pq(j,i,igcm_n2)*mmean(j,i)/mmol(igcm_n2)
339               endif
340
341               !Mixing ratio to density
342               co2(i)=co2(i)*nt                ! CO2 density in cm-3
343               o3p(i)=o3p(i)*nt                ! O3p density in cm-3
344               n2co(i)=n2co(i)*nt              ! N2+CO in cm-3
345c molecular populations
346               n1 = co2(i) * imr1
347               n2 = co2(i) * imr2
348               co2t = n1 + n2
349
350c intermediate collisional rates
351               tt = pt(j,i)*pt(j,i)
352
353               if (pt(j,i).le.175.0) then
354                  k19xca = 3.3e-15
355                  k19xcb = 7.6e-16
356               else
357                  k19xca = 4.2e-12 * exp( -2988.0/pt(j,i) + 303930.0/tt)
358                  k19xcb = 2.1e-12 * exp( -2659.0/pt(j,i) + 223052.0/tt)
359               endif
360               k19xca = k19xca * rfvt
361               k19xcb = k19xcb * rfvt
362               k19cap1 = k19xca * 2.0 * exp( -ee*nu1/pt(j,i) )
363               k19cap2 = k19xca * 2.0 * exp( -ee*nu2/pt(j,i) )
364               k19cbp1 = k19xcb * 2.0 * exp( -ee*nu1/pt(j,i) )
365               k19cbp2 = k19xcb * 2.0 * exp( -ee*nu2/pt(j,i) )
366               d19c = k19xca*co2t + k19xcb*n2co(i)
367               d19cp1 = k19cap1*co2t + k19cbp1*n2co(i)
368               d19cp2 = k19cap2*co2t + k19cbp2*n2co(i)
369                                !
370               k20xc = 3.e-12 * rfvto3p
371               k20cp1 = k20xc * 2.0 * exp( -ee/pt(j,i) * nu1 )
372               k20cp2 = k20xc * 2.0 * exp( -ee/pt(j,i) * nu2 )
373                                !
374               k21xc = 2.49e-11 * 0.5 * rfvv
375               k21cp2 = k21xc * exp( - ee/pt(j,i) * (nu2-nu1) )
376                                !
377               l1 = d19c + k20xc*o3p(i) + k21cp2*n2
378               p1 = ( d19cp1 + k20cp1*o3p(i) ) * n1
379               p12 = k21xc*n1
380                                !
381               l2 = d19c + k20xc*o3p(i) + k21xc*n1
382               p2 = ( d19cp2 + k20cp2*o3p(i) ) * n2
383               p21 = k21cp2*n2
384
385c radiative rates
386               ae1 = 1.3546 * 1.66 / 4.0 * escf1(i)
387               ae2 = ( 1.3452 + 1.1878 ) * 1.66 / 4.0 * escf2(i)
388               l1 = l1 + ae1
389               l2 = l2 + ae2
390
391c solving the system
392               c1 = gamma*nu1**3. * 0.5
393               c2 = gamma*nu2**3. * 0.5
394               a1 = c1 * p1 / (n1*l1)
395               a2 = c2 * p2 / (n2*l2)
396               a12 = (nu1/nu2)**3. * n2/n1 * p12/l1
397               a21 = (nu2/nu1)**3. * n1/n2 * p21/l2
398               el2 = (a2 + a21 * a1 ) / ( 1.0 - a21 * a12 )
399               el1 = a1 + a12 * el2
400               pl1 = el1 * n1 / c1
401               pl2 = el2 * n2 / c2
402
403c  heating rate
404               hr1 = - hplanck*vlight * nu1 * ae1 * pl1
405               hr2 = - hplanck*vlight * nu2 * ae2 * pl2
406               hr(i) = hr1 + hr2
407               dtnlte(j,i)=0.1*hr(i)*pt(j,i)/(4.4*pyy(i)) ! dtnlte in K s-1
408c              write(7,25)pxx(i),hr1,hr2,hr(i),qt
409c  25         format(' ',1p5e12.4)
410
411            endif
412
413         enddo  ! end loop over layers
414      enddo     ! end loop over grid points
415c     close(7)
416c
417
418      end subroutine nltecool
419
420c***********************************************************************
421
422      subroutine interp1(escout,p,nlayer,escin,pin,nl)
423C
424C subroutine to perform linear interpolation in pressure from 1D profile
425C escin(nl) sampled on pressure grid pin(nl) to profile
426C escout(nlayer) on pressure grid p(nlayer).
427C
428      real escout(nlayer),p(nlayer)
429      real escin(nl),pin(nl),wm,wp
430      integer nl,nlayer,n1,n,nm,np
431      do n1=1,nlayer
432         if(p(n1) .gt. 3.5 .or. p(n1) .lt. 4.0e-6) then
433            escout(n1) = 0.0
434         else
435            do n = 1,nl-1
436               if (p(n1).le.pin(n).and.p(n1).ge.pin(n+1)) then
437                  nm=n
438                  np=n+1
439                  wm=abs(pin(np)-p(n1))/(pin(nm)-pin(np))
440                  wp=1.0 - wm
441               endif
442            enddo
443            escout(n1) = escin(nm)*wm + escin(np)*wp
444         endif
445      enddo
446
447      end subroutine interp1
448
449c***********************************************************************
450
451      subroutine interp3(esco1,esco2,esco3,p,nlayer,
452     1     esci1,esci2,esci3,pin,nl)
453C
454C subroutine to perform 3 simultaneous linear interpolations in pressure from
455C 1D profiles esci1-3(nl) sampled on pressure grid pin(nl) to 1D profiles
456C esco1-3(nlayer) on pressure grid p(ngrid,nlayer).
457C
458      real esco1(nlayer),esco2(nlayer),esco3(nlayer),p(nlayer)
459      real esci1(nl),    esci2(nl),    esci3(nl),    pin(nl),wm,wp
460      integer nl,nlayer,n1,n,nm,np
461      do n1=1,nlayer
462         if (p(n1).gt. 3.5 .or. p(n1) .lt. 4.0e-6) then
463            esco1(n1)=0.0
464            esco2(n1)=0.0
465            esco3(n1)=0.0
466         else
467            do n = 1,nl-1
468               if (p(n1).le.pin(n).and.p(n1).ge.pin(n+1)) then
469                  nm=n
470                  np=n+1
471                  wm=abs(pin(np)-p(n1))/(pin(nm)-pin(np))
472                  wp=1.0 - wm
473               endif
474            enddo
475            esco1(n1) = esci1(nm)*wm + esci1(np)*wp
476            esco2(n1) = esci2(nm)*wm + esci2(np)*wp
477            esco3(n1) = esci3(nm)*wm + esci3(np)*wp
478         endif
479      enddo
480
481      end subroutine interp3
482
483      END MODULE nltecool_mod
Note: See TracBrowser for help on using the repository browser.