source: trunk/LMDZ.MARS/libf/phymars/nlte_setup.F @ 1635

Last change on this file since 1635 was 757, checked in by emillour, 12 years ago

Mars GCM:

  • Improvement of the NLTE 15um scheme (for running with nltemodel = 2); now MUCH faster than previously (by a factor 5 or so):
  • Improvements included to the parameterization:
    • Cool-to-space calculation included above P(atm)=1e-10, with a soft merging to the full result (without the CTS approximation) below that level
    • exhaustive cleaning of the code, including FTNCHK and reordering of loops, subroutines and internal calls
    • simplification of the precomputed tables of CO2 bands' atmospheric transmittances
    • the two internal grids (the one used in the CTS region and the one below) have been , in order to reduce the CPU time consumption
    • reading of the spectroscopic histograms is made only once, at the beginning of the GCM, to avoid repetitive readings of ASCII files
    • F90 matrix operations (matmul,...) included.
  • Changes in routines:
    • removed nlte_leedat.F
    • updated nlte_calc.F, nlte_tcool.F, nlte_aux.F
    • updated nlte_commons.h, nlte_paramdef.h
    • added nlte_setup.F
  • Important: The input files (in the NLTEDAT directory) read as input by these routines have changed. the NLTEDAT directory should now on contain only the following files:

deltanu26.dat enelow27.dat hid26-3.dat parametp_Tstar_IAA1204.dat
deltanu27.dat enelow28.dat hid26-4.dat parametp_VC_IAA1204.dat
deltanu28.dat enelow36.dat hid27-1.dat
deltanu36.dat hid26-1.dat hid28-1.dat
enelow26.dat hid26-2.dat hid36-1.dat

FGG+MALV
File size: 6.7 KB
Line 
1c***********************************************************************
2
3      subroutine nlte_setup
4
5c     malv    Oct 09          Adapt mz1d_onlyTCR_MUCHASveces.f to "V09"
6c     malv    Sep 07          Add LU deccomp & repetition option to test CPU
7c     malv    Jan 07          Add new vertical fine-grid for NLTE
8c     apr 06  malv            Read date,effuv from Driver. T fixed at zbott.
9c     2003    fgg             Double precission in UV, Photoq, Conduct & Diff
10c     oct 02  malv            V02: New scheme to allow for continuity eq.
11c     dec 01  malv            See changes/progress of the code in mz1d.actual
12c     nov 01  malv            adapt for parameterizations of tcr y shr
13c     nov 98  malv            add chemical & photochem. processes
14c     jul 98  malv            transic hiperb con zs fuera de la region
15c     equil hidrostatico. smoothing en cr y sh
16c     jan 98    malv            first version
17c***********************************************************************
18
19      implicit none
20
21      include   'datafile.h'
22      include   'nlte_paramdef.h'
23      include   'nlte_commons.h'
24
25
26c***************
27
28c     local variables
29
30      integer   i, k, lun1, lun2
31      real*8          xx
32      character isotcode*2
33
34
35c     formats
36 132  format (i2)
37
38c**********************************************************************
39
40c     *** Groups old 1-d model subroutines SETTINGS and LeeESCTVCISO_dlvr11
41c     *** Both were called in old NLTEdlvr11_SETUP ***
42
43c     *** Old SETTINGS ***
44
45      lun1 = 1
46      lun2 = 2
47
48      do k=1,nisot
49         write (isotcode,132) indexisot(k)
50         open (lun1,
51     $        file=trim(datafile)//'/NLTEDAT/enelow'
52     $        //isotcode//'.dat',status='old')
53         open (lun2,
54     $        file=trim(datafile)//'/NLTEDAT/deltanu'
55     $        //isotcode//'.dat',status='old')
56         read (lun1,*)
57         read (lun2,*)
58         read (lun1,*) (elow(k,i), i=1,nb)
59         read (lun2,*) (deltanu(k,i), i=1,nb)
60         close (lun1)
61         close (lun2)
62      end do
63
64      a1_010_000 = 1.3546d00
65      a2_010_000 = 1.3452d00
66      a3_010_000 = 1.1878d00
67      a4_010_000 = 1.2455d00
68      a1_020_010 = 4.35d0
69     
70
71c     *** Old LeeESCTVCISO_dlvr11 ***
72
73      open( 11, file=trim(datafile)//
74     $     '/NLTEDAT/parametp_Tstar_IAA1204.dat' )
75      read (11, *)
76      do i=1,nztabul
77         read (11,*) lnpnbtab(i), tstar11tab(i),
78     $        tstar21tab(i), tstar31tab(i), tstar41tab(i)
79      enddo
80      close (11)
81
82      open( 12, file=trim(datafile)//
83     $     '/NLTEDAT/parametp_VC_IAA1204.dat' )
84      read (12, *)
85      do i=1,nztabul
86         read (12,*) xx, vc210tab(i), vc310tab(i), vc410tab(i)
87      enddo
88      close (12)
89      xx=xx
90
91
92
93      call LeeHISTOGRMS
94
95c     end subroutine
96
97      return
98      end
99
100
101
102c***********************************************************************
103      subroutine LeeHISTOGRMS
104c***********************************************************************
105
106      implicit none
107
108      include 'datafile.h'
109      include   'nlte_paramdef.h'
110      include   'nlte_commons.h'
111
112
113c     local variables and constants
114      integer   ihist
115
116
117c***********************************************************************
118
119                                ! Banda fundamental
120                                !
121      hisfile = trim(datafile)//
122     $     '/NLTEDAT/hid26-1.dat'
123      ihist = 1
124      call rhist_03 (ihist)
125
126
127                                ! First Hot bands       
128                                !     
129      hisfile = trim(datafile)//
130     $     '/NLTEDAT/hid26-2.dat'
131      ihist = 2
132      call rhist_03 (ihist)
133
134      hisfile = trim(datafile)//
135     $     '/NLTEDAT/hid26-3.dat'
136      ihist = 3
137      call rhist_03 (ihist)
138
139      hisfile = trim(datafile)//
140     $     '/NLTEDAT/hid26-4.dat'
141      ihist = 4
142      call rhist_03 (ihist)
143
144
145
146
147      return
148      end
149
150
151c     *** Old GETK_dlvr11.f ***
152
153c***********************************************************************
154
155      subroutine GETK_dlvr11 (tt)
156
157c***********************************************************************
158
159      implicit none
160
161      include 'nlte_paramdef.h'
162      include 'nlte_commons.h'
163
164c     arguments
165      real              tt      ! i. temperature
166
167!     ! local variables:
168      real*8 k20x, k20xb, k20xc
169      real*8 k19xca,k19xcb,k19xcc
170      real*8 k19xba,k19xbb,k19xbc
171      real*8 k21x,k21xa,k21xb,k21xc
172      real*8 anu, factor , tdt
173      integer   i
174
175c***********************************************************************
176
177      tdt = dble(tt)
178
179                                !! k19 & k20
180
181      k20x = 3.d-12
182      k20xc = k20x * rf20
183      k20xb = 2.d0 * k20xc
184
185      k19xca = 4.2d-12 * exp( -2988.d0/tdt + 303930.d0/tdt**2.d0 )
186      if (tt.le.175.) k19xca = 3.3d-15
187      k19xcb = 2.1d-12 * exp( -2659.d0/tdt + 223052.d0/tdt**2.d0 )
188      if (tt.le.175.) k19xcb = 7.6d-16
189      k19xca = k19xca * rf19
190      k19xcb = k19xcb * rf19
191      k19xcc = k19xcb
192
193      factor = 2.5d0
194      k19xba = factor * k19xca
195      k19xbb = factor * k19xcb
196      k19xbc = factor * k19xcc
197
198      do i = 1, nisot
199
200         k19ba(i) = k19xba
201         k19ca(i) = k19xca
202         k19bb(i) = k19xbb
203         k19cb(i) = k19xcb
204         k19bc(i) = k19xbc
205         k19cc(i) = k19xcc
206
207         k20b(i) = k20xb
208         k20c(i) = k20xc
209
210         anu = dble( nu(i,2)-nu(i,1) )
211
212         k19bap(i) = k19ba(i) * 2.d0 * exp( -ee*anu/tdt )
213         k19bbp(i) = k19bb(i) * 2.d0 * exp( -ee*anu/tdt )
214         k19bcp(i) = k19bc(i) * 2.d0 * exp( -ee*anu/tdt )
215
216         k20bp(i) = k20b(i)*4.d0/2.d0 * exp( -ee/tdt * anu )
217
218         anu = dble( nu(i,1) )
219
220         k19cap(i) = k19ca(i) * 2.d0 * exp( -ee*anu/tdt )
221         k19cbp(i) = k19cb(i) * 2.d0 * exp( -ee*anu/tdt )
222         k19ccp(i) = k19cc(i) * 2.d0 * exp( -ee*anu/tdt )
223
224         k20cp(i) = k20c(i)*2.d0/1.d0 * exp( -ee/tdt * anu )
225
226      end do
227
228
229                                !! k21 &  k23k21c &  k24k21c & k34k21c
230
231      k21x = 2.49d-11
232      k21xb = k21x
233      k21xa = 3.d0/2.d0 * k21xb
234      k21xc = k21xb / 2.d0
235
236      k21xa = k21xa * rf21a
237      k21xb = k21xb * rf21b
238      k21xc = k21xc * rf21c
239
240      do i = 1, nisot
241         k21b(i) = k21xb
242         k21c(i) = k21xc
243         k21bp(i) = k21b(i) *
244     @        exp( -ee/tdt* dble(nu(i,2)-nu(i,1)-nu(1,1)) )
245         k21cp(i) = k21c(i) *
246     @        exp( -ee/tdt * dble(nu(i,1)-nu(1,1)) )
247      end do
248
249      k23k21c = k21xc
250      k24k21c = k21xc
251      k34k21c = k21xc
252      k23k21cp = k23k21c*2.d0/2.d0 *
253     @     exp( -ee/tdt* dble(nu(2,1)-nu(3,1)) )
254      k24k21cp = k24k21c*2.d0/2.d0 *
255     @     exp( -ee/tdt* dble(nu(2,1)-nu(4,1)) )
256      k34k21cp = k34k21c*2.d0/2.d0 *
257     @     exp( -ee/tdt* dble(nu(3,1)-nu(4,1)) )
258
259
260                                !! k33
261
262      k33c = k21x * rf33bc
263      do i=2,nisot
264         k33cp(i) = k33c *
265     @        exp( -ee/tdt * dble(nu(1,2)-nu(1,1)-nu(i,1)) )
266      end do
267
268
269      return
270      end
271
272
273
274
275
276
277
Note: See TracBrowser for help on using the repository browser.