source: trunk/LMDZ.MARS/libf/phymars/mztf_correccion.F @ 455

Last change on this file since 455 was 414, checked in by aslmd, 14 years ago

LMDZ.MARS : NEW NLTE MODEL FROM GRANADA AMIGOS

23/11/11 == FGG + MALV

New parameterization of the NLTE 15 micron cooling. The old parameterization is kept as an option, including or not variable atomic oxygen concentration. A new flag is introduced in callphys.def, nltemodel, to select which parameterization wants to be used (new one, old one with variable [O], or old one with fixed [O], see below). Includes many new subroutines and commons in phymars. Some existing routines are also modified:

-physiq.F. Call to the new subroutine NLTE_leedat in first call. Call to nltecool modified to allow for variable atomic oxygen. Depending on the value of nltemodel, the new subroutine NLTEdlvr09_TCOOL is called instead of nltecool.

-inifis.F. Reading of nltemodel is added.

-callkeys.h Declaration of nltemodel is added.

The following lines should be added to callphys.def (ideally after setting callnlte):

# NLTE 15um scheme to use.
# 0-> Old scheme, static oxygen
# 1-> Old scheme, dynamic oxygen
# 2-> New scheme
nltemodel = 2

A new directory, NLTEDAT, has to be included in datagcm.

Improvements into NLTE NIR heating parameterization to take into account variability with O/CO2 ratio and SZA. A new subroutine, NIR_leedat.F, and a new common, NIRdata.h, are included. A new flag, nircorr, is added in callphys.def, to include or not these corrections. The following files are modified:

-nirco2abs.F: nq and pq are added in the arguments. The corrections factors are interpolated to the GCM grid and included in the clculation. A new subroutine, interpnir, is added at the end of the file.

-physiq.F: Call to NIR_leedat added at first call. Modified call to nirco2abs

-inifis: Reading new flag nircorr.

-callkeys.h: Declaration of nircorr.

The following lines have to be added to callphys.def (ideally after callnirco2):

# NIR NLTE correction for variable SZA and O/CO2?
# matters only if callnirco2=T
# 0-> no correction
# 1-> include correction
nircorr=1

A new data file, NIRcorrection_feb2011.dat, has to be included in datagcm.

Small changes to the molecular diffusion scheme to fix the number of species considered, to avoid problems when compiling with more than 15 tracers (for example, when CH4 is included). Modified subroutines: aeronomars/moldiff.F and aeronomars/moldiffcoeff.F
File size: 15.5 KB
Line 
1c***********************************************************************
2                                               
3        subroutine mztf_correccion (coninf, con, ib, isot, icurt_pop) 
4                                               
5c including the dependence of the absort. coeff. on temp., vibr. temp.,
6c function, etc.., when neccessary. imr is already corrected in his.for
7c we follow pg.39b-43a (l5):                   
8c  tvt1 is the vibr temp of the upper level     
9c  tvt  is the vibr temp of the transition itself           
10c  tvtbs is the vibr temp of the bending mode (used in qv) 
11c  for fundamental bands, they are not used at the moment. 
12c  for the 15 fh and sh bands, only tvt0 is used at the moment.         
13c  for the laser band, all of them are used following pg. 41a -l5- :   
14c    we need s(z) and we can read s(tk) from the histogram (also called
15c    what we have to calculate now is the factor s(z)/s(tk) or following
16c    l5 notebook notation, s_nlte/s_lte.       
17c           s_nlte/s_lte = xfactor = xlower * xqv * xes     
18                                               
19c  icurt_pop = 30 -> Output of populations of the 0200,0220,1000 states
20c            = otro -> no output of these populations
21
22c       oct 92          malv                   
23c       jan 98          malv            version for mz1d         
24c       jul 2011        malv+fgg        adapted to LMD-MGCM
25c***********************************************************************
26                                               
27        implicit none                                 
28                                               
29        include 'nltedefs.h'   
30        include 'nlte_atm.h'
31        include 'nlte_data.h'
32        include 'nlte_results.h'   
33        include 'nlte_curtis.h'       
34                                               
35c arguments                                     
36        integer         ib, isot                             
37        integer         icurt_pop           ! output of Fermi states population
38        real*8          con(nzy), coninf                       
39                                               
40! local variables                               
41        integer         i                                     
42        real*8  tvt0(nzy),tvt1(nzy),tvtbs(nzy), zld(nl),zyd(nzy)               
43        real    xalfa, xbeta, xtv1000, xtv0200, xtv0220, xfactor     
44        real    xqv, xnu_trans, xtv_trans, xes, xlower   
45c***********************************************************************
46                                 
47        xfactor = 1.0
48
49        do i=1,nzy
50          zyd(i) = dble(zy(i))
51        enddo
52        do i=1,nl                                     
53          zld(i) = dble( zl(i) )               
54        end do                                 
55                                               
56! tvtbs is the bending mode of the molecule. used in xqv.   
57        if (isot.eq.1) call interdp (tvtbs,zyd,nzy, v626t1,zld,nl, 1 ) 
58        if (isot.eq.2) call interdp (tvtbs,zyd,nzy, v628t1,zld,nl, 1 ) 
59        if (isot.eq.3) call interdp (tvtbs,zyd,nzy, v636t1,zld,nl, 1 ) 
60        if (isot.eq.4) call interdp (tvtbs,zyd,nzy, v627t1,zld,nl, 1 ) 
61        if (isot.eq.5) call interdp (tvtbs,zyd,nzy, vcot1,zld,nl, 1 )   
62
63! tvt0 is the lower level of the transition. used in xlower.           
64        if (ib.eq.2 .or. ib.eq.3 .or. ib.eq.4 .or. ib.eq.15) then 
65          if (isot.eq.1) call interdp (tvt0,zyd,nzy, v626t1,zld,nl, 1 )
66          if (isot.eq.2) call interdp (tvt0,zyd,nzy, v628t1,zld,nl, 1 )
67          if (isot.eq.3) call interdp (tvt0,zyd,nzy, v636t1,zld,nl, 1 )
68          if (isot.eq.4) call interdp (tvt0,zyd,nzy, v627t1,zld,nl, 1 )
69        elseif (ib.eq.6 .or. ib.eq.8 .or. ib.eq.10     
70     @          .or. ib.eq.13 .or. ib.eq.14                 
71     @          .or. ib.eq.17 .or. ib.eq.19 .or. ib.eq.20) then         
72          if (isot.eq.1) call interdp ( tvt0,zyd,nzy, v626t2,zld,nl, 1 )
73          if (isot.eq.2) call interdp ( tvt0,zyd,nzy, v628t2,zld,nl, 1 )
74          if (isot.eq.3) call interdp ( tvt0,zyd,nzy, v636t2,zld,nl, 1 )
75          if (isot.eq.4) then 
76                call interdp ( tvt0,zyd,nzy, v627t2,zld,nl, 1 )
77          endif                                       
78        else                                           
79          do i=1,nzy                                   
80            tvt0(i) = dble( ty(i) )                   
81          end do                                       
82        end if                                         
83                                               
84c tvt is the vt of the transition. used in xes.
85c since xes=1.0 except for the laser bands, tvt is only needed for  them
86c but it is actually calculated from the tv of the upper and lower level
87c of the transition. hence, only tvt1 remains to be read for the laser b
88c tvt1 is the upper level of the transition.   
89        if (ib.eq.13 .or. ib.eq.14) then
90          if (isot.eq.1) call interdp ( tvt1,zyd,nzy, v626t4,zld,nl, 1 )
91          if (isot.eq.2) call interdp ( tvt1,zyd,nzy, v628t4,zld,nl, 1 )
92          if (isot.eq.3) call interdp ( tvt1,zyd,nzy, v636t4,zld,nl, 1 )
93          if (isot.eq.4) call interdp ( tvt1,zyd,nzy, v627t4,zld,nl, 1 )
94        end if
95                                               
96c  here we weight the absorber amount by a factor which compensate the l
97c  value of the strength read from hitran. we use that factor in order t
98c  correct the product s*m when we later multiply those two variables. 
99                                               
100!        if ( isot.eq.1 .and. icurt_pop.eq.30 ) then
101!           open (30, file='020populations.dat')
102!           write (30,*) ' z  tv(020) tv0200 tv0220 tv1000 '
103!        endif
104
105        do i=1,nzy                                     
106                                               
107          if (isot.eq.1) then                 
108
109            !!! vt of the 3 levels in (020)  (see pag. 36a-sn1 for this)       
110            xalfa = 1.d0/2.d0* exp( dble(-ee*(nu12_1000-nu(1,2))/ty(i)) )     
111            xbeta = 1.d0/2.d0* exp( dble(-ee*(nu12_0200-nu(1,2))/ty(i)) )     
112            xtv0200 = dble( - ee * nu12_0200 ) /       
113     @        ( log( xbeta/(1.d0+xalfa+xbeta) ) -
114     @                      dble(ee*nu(1,2))/tvt0(i) )   
115            xtv0220 = dble( - ee * nu(1,2) ) /         
116     @        ( log( 1.d0/(1.d0+xalfa+xbeta) ) -
117     @                      dble(ee*nu(1,2))/tvt0(i) )     
118            xtv1000 = dble( - ee * nu12_1000 ) /       
119     @        ( log( xalfa/(1.d0+xalfa+xbeta) ) -
120     @                      dble(ee*nu(1,2))/tvt0(i) )
121            !!! correccion 8-Nov-04 (see pag.9b-Marte4-)
122            xtv0200 = dble( - ee * nu12_0200 /       
123     @       ( log(4.*xbeta/(1.d0+xalfa+xbeta)) - ee*nu(1,2)/tvt0(i) ) )   
124            xtv0220 = dble( - ee * nu(1,2) /         
125     @        ( log(2./(1.d0+xalfa+xbeta)) - ee*nu(1,2)/tvt0(i) ) )     
126            xtv1000 = dble( - ee * nu12_1000 /       
127     @       ( log(4.*xalfa/(1.d0+xalfa+xbeta)) - ee*nu(1,2)/tvt0(i) ) )
128
129!            if ( icurt_pop.eq.30 ) then
130!               write (30,'( 1x,f7.2, 3x,f8.3, 2x,3(1x,f8.3) )')
131!     @           zx(i), tvt0(i), xtv0200, xtv0220, xtv1000
132!            endif
133               
134            !!! xlower and xes for the band           
135            if (ib.eq.19) then                         
136              xlower = exp( dble(ee*elow(isot,ib)) *   
137     @                          ( 1.d0/dble(ty(i))-1.d0/xtv0200 ) )       
138              xes = 1.0d0                             
139            elseif (ib.eq.17) then                     
140              xlower = exp( dble(ee*elow(isot,ib)) *   
141     @                          ( 1.d0/dble(ty(i))-1.d0/xtv1000 ) )       
142              xes = 1.0d0                             
143            elseif (ib.eq.20) then                     
144              xlower = exp( dble(ee*elow(isot,ib)) *   
145     @                          ( 1.d0/dble(ty(i))-1.d0/xtv0220 ) )       
146              xes = 1.0d0                             
147            elseif (ib.eq.14) then                     
148              xlower = exp( dble(ee*nu12_1000) *       
149     @                          ( 1.d0/dble(ty(i))-1.d0/xtv1000 ) )       
150              xnu_trans = dble( nu(1,4)-nu12_1000 )   
151              xtv_trans = xnu_trans / dble(nu(1,4)/tvt1(i)-
152     @                      nu12_1000/xtv1000) 
153              xes = (1.d0-exp( dble(-ee*xnu_trans/xtv_trans) )) / 
154     @                  (1.d0-exp( dble(-ee*xnu_trans/ty(i)) ))     
155            elseif (ib.eq.13) then                     
156              xlower = exp( dble(ee*nu12_0200) *       
157     @                          ( 1.d0/dble(ty(i))-1.d0/xtv0200 ) )       
158              xnu_trans = dble(nu(1,4)-nu12_0200)     
159              xtv_trans = xnu_trans / dble(nu(1,4)/tvt1(i)-
160     @                      nu12_0200/xtv0200) 
161              xes = (1.d0-exp( dble(-ee*xnu_trans/xtv_trans) )) / 
162     @                  (1.d0-exp( dble(-ee*xnu_trans/ty(i)) ))     
163            else                                       
164              xlower = exp( dble(ee*elow(isot,ib)) *   
165     @                          ( 1.d0/dble(ty(i))-1.d0/tvt0(i) ) )       
166              xes = 1.0d0                             
167            end if                                     
168            xqv = (1.d0-exp( dble(-ee*667.3801/tvtbs(i)) )) /     
169     @                (1.d0-exp( dble(-ee*667.3801/ty(i)) ))
170            xfactor = xlower * xqv**2.d0 * xes         
171                                                       
172          elseif (isot.eq.2) then                     
173                                               
174            xalfa = 1.d0/2.d0* exp( dble(-ee*(nu22_1000-nu(2,2))/
175     @                      ty(i)) )     
176            xbeta = 1.d0/2.d0* exp( dble(-ee*(nu22_0200-nu(2,2))/
177     @                      ty(i)) )     
178            xtv0200 = dble( - ee * nu22_0200 ) /       
179     @        ( log( xbeta/(1.d0+xalfa+xbeta) ) - dble(ee*nu(2,2))/
180     @                      tvt0(i) )   
181            xtv1000 = dble( - ee * nu22_1000 ) /       
182     @        ( log( xalfa/(1.d0+xalfa+xbeta) ) - dble(ee*nu(2,2))/
183     @                      tvt0(i) )   
184                                               
185            if (ib.eq.14) then                         
186              xlower = exp( dble(ee*nu22_1000) *       
187     @                          ( 1.d0/dble(ty(i))-1.d0/xtv1000 ) )       
188              xnu_trans = dble(nu(2,4)-nu22_1000)     
189              xtv_trans = xnu_trans / dble(nu(2,4)/tvt1(i)-nu22_1000/
190     @                      xtv1000) 
191              xes = (1.d0-exp( dble(-ee*xnu_trans/xtv_trans) )) / 
192     @                  (1.d0-exp( dble(-ee*xnu_trans/ty(i)) ))     
193            elseif (ib.eq.13) then                     
194              xlower = exp( dble(ee*nu22_0200) *       
195     @                          ( 1.d0/dble(ty(i))-1.d0/xtv0200 ) )       
196              xnu_trans = dble( nu(2,4)-nu22_0200 )   
197              xtv_trans = xnu_trans / dble(nu(2,4)/tvt1(i)-nu22_0200/
198     @                      xtv0200) 
199              xes = (1.d0-exp( dble(-ee*xnu_trans/xtv_trans) )) / 
200     @                  (1.d0-exp( dble(-ee*xnu_trans/ty(i)) ))     
201            else                                       
202              xlower = exp( dble(ee*elow(isot,ib)) *   
203     @                          ( 1.d0/dble(ty(i))-1.d0/tvt0(i) ) )       
204              xes = 1.0d0                             
205            end if                                     
206            xqv = (1.d0-exp( dble(-ee*662.3734/tvtbs(i)) )) /     
207     @                (1.d0-exp( dble(-ee*662.3734/ty(i))  ))           
208            xfactor = xlower * xqv**2.d0 * xes         
209                                               
210          elseif (isot.eq.3) then                     
211                                               
212            xalfa = 1.d0/2.d0* exp( dble(-ee*(nu32_1000-nu(3,2))/
213     @                      ty(i)) )     
214            xbeta = 1.d0/2.d0* exp( dble(-ee*(nu32_0200-nu(3,2))/
215     @                      ty(i)) )     
216            xtv0200 = dble( - ee * nu32_0200 ) /       
217     @        ( log( xbeta/(1.d0+xalfa+xbeta) ) - dble(ee*nu(3,2))/
218     @                      tvt0(i) )   
219            xtv1000 = dble( - ee * nu32_1000 ) /       
220     @        ( log( xalfa/(1.d0+xalfa+xbeta) ) - dble(ee*nu(3,2))/
221     @                      tvt0(i) )   
222                                               
223            if (ib.eq.14) then                         
224              xlower = exp( dble(ee*nu32_1000) *       
225     @                          ( 1.d0/dble(ty(i))-1.d0/xtv1000 ) )       
226              xnu_trans = dble( nu(3,4)-nu32_1000 )   
227              xtv_trans = xnu_trans / dble(nu(3,4)/tvt1(i)-nu32_1000/
228     @                      xtv1000) 
229              xes = (1.d0-exp( dble(-ee*xnu_trans/xtv_trans) )) / 
230     @                  (1.d0-exp( dble(-ee*xnu_trans/ty(i)) ))     
231            elseif (ib.eq.13) then                     
232              xlower = exp( dble(ee*nu32_0200) *       
233     @                          ( 1.d0/dble(ty(i))-1.d0/xtv0200 ) )       
234              xnu_trans = dble( nu(3,4)-nu32_0200 )   
235              xtv_trans = xnu_trans / dble(nu(3,4)/tvt1(i)-nu32_0200/
236     @                      xtv0200) 
237              xes = (1.d0-exp( dble(-ee*xnu_trans/xtv_trans) )) / 
238     @                  (1.d0-exp( dble(-ee*xnu_trans/ty(i)) ))     
239            else                                       
240              xlower = exp( dble(ee*elow(isot,ib)) *   
241     @                          ( 1.d0/dble(ty(i))-1.d0/tvt0(i) ) )       
242              xes = 1.0d0                             
243            end if                                     
244            xqv = (1.d0-exp( dble(-ee*648.4784/tvtbs(i)) )) /     
245     @                (1.d0-exp( dble(-ee*648.4784/ty(i))  ))           
246            xfactor = xlower * xqv**2.d0 * xes         
247                                               
248          elseif (isot.eq.4) then                     
249                                               
250            xalfa = 1.d0/2.d0* exp( dble(-ee*(nu42_1000-nu(4,2))/
251     @                      ty(i)) )     
252            xbeta = 1.d0/2.d0* exp( dble(-ee*(nu42_0200-nu(4,2))/
253     @                      ty(i)) )     
254            xtv0200 = dble( - ee * nu42_0200 ) /       
255     @        ( log( xbeta/(1.d0+xalfa+xbeta) ) - dble(ee*nu(4,2))/
256     @                      tvt0(i) )   
257            xtv1000 = dble( - ee * nu42_1000 ) /       
258     @        ( log( xalfa/(1.d0+xalfa+xbeta) ) - dble(ee*nu(4,2))/
259     @                      tvt0(i) )   
260                                               
261            if (ib.eq.14) then                         
262              xlower = exp( dble(ee*nu42_1000) *       
263     @                          ( 1.d0/dble(ty(i))-1.d0/xtv1000 ) )       
264              xnu_trans = dble( nu(4,4)-nu42_1000 )   
265              xtv_trans = xnu_trans / dble(nu(4,4)/tvt1(i)-nu42_1000/
266     @                      xtv1000) 
267              xes = (1.d0-exp( dble(-ee*xnu_trans/xtv_trans) )) / 
268     @                  (1.d0-exp( dble(-ee*xnu_trans/ty(i)) ))     
269            elseif (ib.eq.13) then                     
270              xlower = exp( dble(ee*nu42_0200) *       
271     $  ( 1.d0/dble(ty(i))-1.d0/xtv0200 ) )     
272              xnu_trans = dble( nu(4,4)-nu42_0200 )   
273              xtv_trans = xnu_trans / dble(nu(4,4)/tvt1(i)-nu42_0200/
274     @                      xtv0200) 
275              xes = (1.d0-exp( dble(-ee*xnu_trans/xtv_trans) )) / 
276     @                  (1.d0-exp( dble(-ee*xnu_trans/ty(i)) ))     
277            else                                       
278              xlower = exp( dble(ee*elow(isot,ib)) *   
279     @                          ( 1.d0/dble(ty(i))-1.d0/tvt0(i) ) )       
280              xes = 1.0d0                             
281            end if                                     
282            xqv = (1.d0-exp( dble(-ee*664.7289/tvtbs(i)) )) /     
283     @                (1.d0-exp( dble(-ee*664.7289/ty(i))  ))           
284            xfactor = xlower * xqv**2.d0 * xes         
285                                               
286          elseif (isot.eq.5 .and. ib.eq.1) then       
287                                                             
288            xlower = 1.d0                             
289            xes = 1.0d0                               
290            xqv = (1.d0-exp( dble(-ee*nuco_10/tvtbs(i)) )) /         
291     @                (1.d0-exp( dble(-ee*nuco_10/ty(i))  ))   
292            xfactor = xlower * xqv * xes         
293                                               
294          end if                                       
295                                               
296          con(i) = con(i) * xfactor                   
297          if (i.eq.nzy) coninf = coninf * xfactor       
298                                               
299        end do                                         
300                   
301!        if ( isot.eq.1 .and. icurt_pop.eq.30 ) then
302!           close (30)
303!        endif
304                           
305        return                                         
306        end                                           
Note: See TracBrowser for help on using the repository browser.