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

Last change on this file since 3393 was 3016, checked in by emillour, 16 months ago

Mars PCM:
Code cleanup: nltedata.h is only included in nltecool.F so turn this data
into module data there. Also convert lteparams.h into module nlteparams_h.F90.
And while at it also turn nlthermeq.F and blendrad.F into modules.
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      implicit none
196
197      include "callkeys.h"
198
199c Input and output variables
200c
201      integer,intent(in) :: ngrid               ! no. of horiz. gridpoints
202      integer,intent(in) :: nlayer              ! no. of atmospheric layers
203      integer,intent(in) :: nq                  ! no. of tracers
204      real,intent(in) :: pplay(ngrid,nlayer)    ! input pressure grid (Pa)
205      real,intent(in) :: pt(ngrid,nlayer)       ! input temperatures (K)
206      real,intent(in) :: pq(ngrid,nlayer,nq)    ! input mmrs (kg/kg_air)
207      real,intent(out) :: dtnlte(ngrid,nlayer)  ! output temp. tendencies (K/s)
208
209c
210c Standard atmosphere variables
211c
212      real        nt                              ! number density [cm-3]
213      real        co2(nlayer)                     !  "   of CO2
214      real        o3p(nlayer)                     !  "   of atomic oxygen
215      real        n2co(nlayer)                    !  "  of N2 + CO
216      real        pyy(nlayer)                     ! auxiliary pressure grid
217
218c
219c Vectors and indexes for the tabulation of escape functions and VMR
220c
221c                 np                              ! # data points in tabulation
222c                 pnb(np)                         ! Pressure in tabulation
223c                 ef1(np)                         ! Esc.funct.#1, tabulated
224c                 ef2(np)                         ! Esc.funct.#2, tabulated
225c                 co2vmr(np)                      ! CO2 VMR tabulated
226c                 o3pvmr(np)                      ! CO2 VMR tabulated
227c                 n2covmr(np)                     ! N2+CO VMR tabulated
228      real        escf1(nlayer)                   ! Esc.funct.#1, interpolated
229      real        escf2(nlayer)                   ! Esc.funct.#2, interpolated
230
231
232c
233c Local Constants
234c
235      real       nu1, nu2                         ! freq. of energy levels
236      real       imr1, imr2                       ! isotopic abundances
237      real       hplanck, gamma, vlight           ! physical constants
238      real       ee
239      real       rfvt                             ! collisional rate
240      real       rfvto3p                          !     "
241      real       rfvv                             !     "
242
243c
244c Local variables for the main loop
245c
246      real       n1, n2, co2t                     ! ground populations
247      real       l1, p1, p12                      ! prod & losses
248      real       l2, p2, p21
249      real       tt                               ! dummy variable
250      real       c1, c2                           ! molecular constants
251      real       ae1, ae2                         ! einstein spontaneous emission
252      real       a1, a2, a12, a21
253      real       pl1, pl2
254      real       el1, el2
255      real       hr1, hr2                         ! heating rate due to each band
256      real       hr(nlayer)                       ! total heating rate
257
258c
259c Indexes
260c
261      integer    i
262      integer    j,ii
263
264c
265c Rate coefficients
266c
267      real       k19xca, k19xcb
268      real       k19cap1, k19cap2
269      real       k19cbp1, k19cbp2
270      real       d19c, d19cp1, d19cp2
271      real       k20xc, k20cp1, k20cp2
272      real       k21xc, k21cp2
273
274      logical,save :: firstcall=.true.
275!$OMP THREADPRIVATE(firstcall)
276
277c
278c Data
279c
280      data       nu1, nu2, hplanck, gamma, vlight, ee/
281     1     667.38, 662.3734, 6.6261e-27, 1.191e-5, 3.e10, 1.438769/
282
283c*************************************************************************
284c       PROGRAM  STARTS
285c*************************************************************************
286
287      imr1 = 0.987
288      imr2 = 0.00408 + 0.0112
289      rfvt = 0.1
290      rfvto3p = 1.0
291      rfvv = 0.1
292
293      if(firstcall) then
294
295         do i=1,np
296            pnb(i)=1.0e-4*exp(pnb(i)) ! p into Pa
297         end do
298
299         firstcall = .false.
300
301      endif
302c
303c MAIN LOOP, for each gridpoint and altitude:
304c
305      do j=1,ngrid  ! loop over grid points
306c
307c set up local pressure grid
308c
309         do ii=1,nlayer
310            pyy(ii)=pplay(j,ii)
311         enddo
312!
313! Interpolate escape functions and VMR to the desired grid
314!
315         call interp1(escf2,pyy,nlayer,ef2,pnb,np)
316         call interp1(escf1,pyy,nlayer,ef1,pnb,np)
317         if(nltemodel.eq.0) then
318            call interp3(co2,o3p,n2co,pyy,nlayer,
319     &           co2vmr,o3pvmr,n2covmr,pnb,np)
320         endif
321       
322         do i=1,nlayer  ! loop over layers
323C
324C test if p lies outside range (p > 3.5 Pa)
325C changed to 1 Pa since transition will always be higher than this
326C
327            if(pyy(i) .gt. 1.0 .or. pyy(i) .lt. 4.0e-6) then
328               hr(i)=0.0
329               dtnlte(j,i)=0.0
330            else
331c
332c           if(pt(j,i).lt.1.0)print*,pt(j,i)
333               nt = pyy(i)/(1.381e-17*pt(j,i)) ! nt in cm-3
334               !Dynamic composition
335               if(nltemodel.eq.1) then
336                  co2(i)=pq(j,i,igcm_co2)*mmean(j,i)/mmol(igcm_co2)
337                  o3p(i)=pq(j,i,igcm_o)*mmean(j,i)/mmol(igcm_o)
338                  n2co(i)=pq(j,i,igcm_co)*mmean(j,i)/mmol(igcm_co) +
339     $                 pq(j,i,igcm_n2)*mmean(j,i)/mmol(igcm_n2)
340               endif
341
342               !Mixing ratio to density
343               co2(i)=co2(i)*nt                ! CO2 density in cm-3
344               o3p(i)=o3p(i)*nt                ! O3p density in cm-3
345               n2co(i)=n2co(i)*nt              ! N2+CO in cm-3
346c molecular populations
347               n1 = co2(i) * imr1
348               n2 = co2(i) * imr2
349               co2t = n1 + n2
350
351c intermediate collisional rates
352               tt = pt(j,i)*pt(j,i)
353
354               if (pt(j,i).le.175.0) then
355                  k19xca = 3.3e-15
356                  k19xcb = 7.6e-16
357               else
358                  k19xca = 4.2e-12 * exp( -2988.0/pt(j,i) + 303930.0/tt)
359                  k19xcb = 2.1e-12 * exp( -2659.0/pt(j,i) + 223052.0/tt)
360               endif
361               k19xca = k19xca * rfvt
362               k19xcb = k19xcb * rfvt
363               k19cap1 = k19xca * 2.0 * exp( -ee*nu1/pt(j,i) )
364               k19cap2 = k19xca * 2.0 * exp( -ee*nu2/pt(j,i) )
365               k19cbp1 = k19xcb * 2.0 * exp( -ee*nu1/pt(j,i) )
366               k19cbp2 = k19xcb * 2.0 * exp( -ee*nu2/pt(j,i) )
367               d19c = k19xca*co2t + k19xcb*n2co(i)
368               d19cp1 = k19cap1*co2t + k19cbp1*n2co(i)
369               d19cp2 = k19cap2*co2t + k19cbp2*n2co(i)
370                                !
371               k20xc = 3.e-12 * rfvto3p
372               k20cp1 = k20xc * 2.0 * exp( -ee/pt(j,i) * nu1 )
373               k20cp2 = k20xc * 2.0 * exp( -ee/pt(j,i) * nu2 )
374                                !
375               k21xc = 2.49e-11 * 0.5 * rfvv
376               k21cp2 = k21xc * exp( - ee/pt(j,i) * (nu2-nu1) )
377                                !
378               l1 = d19c + k20xc*o3p(i) + k21cp2*n2
379               p1 = ( d19cp1 + k20cp1*o3p(i) ) * n1
380               p12 = k21xc*n1
381                                !
382               l2 = d19c + k20xc*o3p(i) + k21xc*n1
383               p2 = ( d19cp2 + k20cp2*o3p(i) ) * n2
384               p21 = k21cp2*n2
385
386c radiative rates
387               ae1 = 1.3546 * 1.66 / 4.0 * escf1(i)
388               ae2 = ( 1.3452 + 1.1878 ) * 1.66 / 4.0 * escf2(i)
389               l1 = l1 + ae1
390               l2 = l2 + ae2
391
392c solving the system
393               c1 = gamma*nu1**3. * 0.5
394               c2 = gamma*nu2**3. * 0.5
395               a1 = c1 * p1 / (n1*l1)
396               a2 = c2 * p2 / (n2*l2)
397               a12 = (nu1/nu2)**3. * n2/n1 * p12/l1
398               a21 = (nu2/nu1)**3. * n1/n2 * p21/l2
399               el2 = (a2 + a21 * a1 ) / ( 1.0 - a21 * a12 )
400               el1 = a1 + a12 * el2
401               pl1 = el1 * n1 / c1
402               pl2 = el2 * n2 / c2
403
404c  heating rate
405               hr1 = - hplanck*vlight * nu1 * ae1 * pl1
406               hr2 = - hplanck*vlight * nu2 * ae2 * pl2
407               hr(i) = hr1 + hr2
408               dtnlte(j,i)=0.1*hr(i)*pt(j,i)/(4.4*pyy(i)) ! dtnlte in K s-1
409c              write(7,25)pxx(i),hr1,hr2,hr(i),qt
410c  25         format(' ',1p5e12.4)
411
412            endif
413
414         enddo  ! end loop over layers
415      enddo     ! end loop over grid points
416c     close(7)
417c
418
419      end subroutine nltecool
420
421c***********************************************************************
422
423      subroutine interp1(escout,p,nlayer,escin,pin,nl)
424C
425C subroutine to perform linear interpolation in pressure from 1D profile
426C escin(nl) sampled on pressure grid pin(nl) to profile
427C escout(nlayer) on pressure grid p(nlayer).
428C
429      real escout(nlayer),p(nlayer)
430      real escin(nl),pin(nl),wm,wp
431      integer nl,nlayer,n1,n,nm,np
432      do n1=1,nlayer
433         if(p(n1) .gt. 3.5 .or. p(n1) .lt. 4.0e-6) then
434            escout(n1) = 0.0
435         else
436            do n = 1,nl-1
437               if (p(n1).le.pin(n).and.p(n1).ge.pin(n+1)) then
438                  nm=n
439                  np=n+1
440                  wm=abs(pin(np)-p(n1))/(pin(nm)-pin(np))
441                  wp=1.0 - wm
442               endif
443            enddo
444            escout(n1) = escin(nm)*wm + escin(np)*wp
445         endif
446      enddo
447
448      end subroutine interp1
449
450c***********************************************************************
451
452      subroutine interp3(esco1,esco2,esco3,p,nlayer,
453     1     esci1,esci2,esci3,pin,nl)
454C
455C subroutine to perform 3 simultaneous linear interpolations in pressure from
456C 1D profiles esci1-3(nl) sampled on pressure grid pin(nl) to 1D profiles
457C esco1-3(nlayer) on pressure grid p(ngrid,nlayer).
458C
459      real esco1(nlayer),esco2(nlayer),esco3(nlayer),p(nlayer)
460      real esci1(nl),    esci2(nl),    esci3(nl),    pin(nl),wm,wp
461      integer nl,nlayer,n1,n,nm,np
462      do n1=1,nlayer
463         if (p(n1).gt. 3.5 .or. p(n1) .lt. 4.0e-6) then
464            esco1(n1)=0.0
465            esco2(n1)=0.0
466            esco3(n1)=0.0
467         else
468            do n = 1,nl-1
469               if (p(n1).le.pin(n).and.p(n1).ge.pin(n+1)) then
470                  nm=n
471                  np=n+1
472                  wm=abs(pin(np)-p(n1))/(pin(nm)-pin(np))
473                  wp=1.0 - wm
474               endif
475            enddo
476            esco1(n1) = esci1(nm)*wm + esci1(np)*wp
477            esco2(n1) = esci2(nm)*wm + esci2(np)*wp
478            esco3(n1) = esci3(nm)*wm + esci3(np)*wp
479         endif
480      enddo
481
482      end subroutine interp3
483
484      END MODULE nltecool_mod
Note: See TracBrowser for help on using the repository browser.