c*********************************************************************** subroutine mztf_correccion (coninf, con, ib, isot, icurt_pop) c including the dependence of the absort. coeff. on temp., vibr. temp., c function, etc.., when neccessary. imr is already corrected in his.for c we follow pg.39b-43a (l5): c tvt1 is the vibr temp of the upper level c tvt is the vibr temp of the transition itself c tvtbs is the vibr temp of the bending mode (used in qv) c for fundamental bands, they are not used at the moment. c for the 15 fh and sh bands, only tvt0 is used at the moment. c for the laser band, all of them are used following pg. 41a -l5- : c we need s(z) and we can read s(tk) from the histogram (also called c what we have to calculate now is the factor s(z)/s(tk) or following c l5 notebook notation, s_nlte/s_lte. c s_nlte/s_lte = xfactor = xlower * xqv * xes c icurt_pop = 30 -> Output of populations of the 0200,0220,1000 states c = otro -> no output of these populations c oct 92 malv c jan 98 malv version for mz1d c jul 2011 malv+fgg adapted to LMD-MGCM c*********************************************************************** implicit none include 'nltedefs.h' include 'nlte_atm.h' include 'nlte_data.h' include 'nlte_results.h' include 'nlte_curtis.h' c arguments integer ib, isot integer icurt_pop ! output of Fermi states population real*8 con(nzy), coninf ! local variables integer i real*8 tvt0(nzy),tvt1(nzy),tvtbs(nzy), zld(nl),zyd(nzy) real xalfa, xbeta, xtv1000, xtv0200, xtv0220, xfactor real xqv, xnu_trans, xtv_trans, xes, xlower c*********************************************************************** xfactor = 1.0 do i=1,nzy zyd(i) = dble(zy(i)) enddo do i=1,nl zld(i) = dble( zl(i) ) end do ! tvtbs is the bending mode of the molecule. used in xqv. if (isot.eq.1) call interdp (tvtbs,zyd,nzy, v626t1,zld,nl, 1 ) if (isot.eq.2) call interdp (tvtbs,zyd,nzy, v628t1,zld,nl, 1 ) if (isot.eq.3) call interdp (tvtbs,zyd,nzy, v636t1,zld,nl, 1 ) if (isot.eq.4) call interdp (tvtbs,zyd,nzy, v627t1,zld,nl, 1 ) if (isot.eq.5) call interdp (tvtbs,zyd,nzy, vcot1,zld,nl, 1 ) ! tvt0 is the lower level of the transition. used in xlower. if (ib.eq.2 .or. ib.eq.3 .or. ib.eq.4 .or. ib.eq.15) then if (isot.eq.1) call interdp (tvt0,zyd,nzy, v626t1,zld,nl, 1 ) if (isot.eq.2) call interdp (tvt0,zyd,nzy, v628t1,zld,nl, 1 ) if (isot.eq.3) call interdp (tvt0,zyd,nzy, v636t1,zld,nl, 1 ) if (isot.eq.4) call interdp (tvt0,zyd,nzy, v627t1,zld,nl, 1 ) elseif (ib.eq.6 .or. ib.eq.8 .or. ib.eq.10 @ .or. ib.eq.13 .or. ib.eq.14 @ .or. ib.eq.17 .or. ib.eq.19 .or. ib.eq.20) then if (isot.eq.1) call interdp ( tvt0,zyd,nzy, v626t2,zld,nl, 1 ) if (isot.eq.2) call interdp ( tvt0,zyd,nzy, v628t2,zld,nl, 1 ) if (isot.eq.3) call interdp ( tvt0,zyd,nzy, v636t2,zld,nl, 1 ) if (isot.eq.4) then call interdp ( tvt0,zyd,nzy, v627t2,zld,nl, 1 ) endif else do i=1,nzy tvt0(i) = dble( ty(i) ) end do end if c tvt is the vt of the transition. used in xes. c since xes=1.0 except for the laser bands, tvt is only needed for them c but it is actually calculated from the tv of the upper and lower level c of the transition. hence, only tvt1 remains to be read for the laser b c tvt1 is the upper level of the transition. if (ib.eq.13 .or. ib.eq.14) then if (isot.eq.1) call interdp ( tvt1,zyd,nzy, v626t4,zld,nl, 1 ) if (isot.eq.2) call interdp ( tvt1,zyd,nzy, v628t4,zld,nl, 1 ) if (isot.eq.3) call interdp ( tvt1,zyd,nzy, v636t4,zld,nl, 1 ) if (isot.eq.4) call interdp ( tvt1,zyd,nzy, v627t4,zld,nl, 1 ) end if c here we weight the absorber amount by a factor which compensate the l c value of the strength read from hitran. we use that factor in order t c correct the product s*m when we later multiply those two variables. ! if ( isot.eq.1 .and. icurt_pop.eq.30 ) then ! open (30, file='020populations.dat') ! write (30,*) ' z tv(020) tv0200 tv0220 tv1000 ' ! endif do i=1,nzy if (isot.eq.1) then !!! vt of the 3 levels in (020) (see pag. 36a-sn1 for this) xalfa = 1.d0/2.d0* exp( dble(-ee*(nu12_1000-nu(1,2))/ty(i)) ) xbeta = 1.d0/2.d0* exp( dble(-ee*(nu12_0200-nu(1,2))/ty(i)) ) xtv0200 = dble( - ee * nu12_0200 ) / @ ( log( xbeta/(1.d0+xalfa+xbeta) ) - @ dble(ee*nu(1,2))/tvt0(i) ) xtv0220 = dble( - ee * nu(1,2) ) / @ ( log( 1.d0/(1.d0+xalfa+xbeta) ) - @ dble(ee*nu(1,2))/tvt0(i) ) xtv1000 = dble( - ee * nu12_1000 ) / @ ( log( xalfa/(1.d0+xalfa+xbeta) ) - @ dble(ee*nu(1,2))/tvt0(i) ) !!! correccion 8-Nov-04 (see pag.9b-Marte4-) xtv0200 = dble( - ee * nu12_0200 / @ ( log(4.*xbeta/(1.d0+xalfa+xbeta)) - ee*nu(1,2)/tvt0(i) ) ) xtv0220 = dble( - ee * nu(1,2) / @ ( log(2./(1.d0+xalfa+xbeta)) - ee*nu(1,2)/tvt0(i) ) ) xtv1000 = dble( - ee * nu12_1000 / @ ( log(4.*xalfa/(1.d0+xalfa+xbeta)) - ee*nu(1,2)/tvt0(i) ) ) ! if ( icurt_pop.eq.30 ) then ! write (30,'( 1x,f7.2, 3x,f8.3, 2x,3(1x,f8.3) )') ! @ zx(i), tvt0(i), xtv0200, xtv0220, xtv1000 ! endif !!! xlower and xes for the band if (ib.eq.19) then xlower = exp( dble(ee*elow(isot,ib)) * @ ( 1.d0/dble(ty(i))-1.d0/xtv0200 ) ) xes = 1.0d0 elseif (ib.eq.17) then xlower = exp( dble(ee*elow(isot,ib)) * @ ( 1.d0/dble(ty(i))-1.d0/xtv1000 ) ) xes = 1.0d0 elseif (ib.eq.20) then xlower = exp( dble(ee*elow(isot,ib)) * @ ( 1.d0/dble(ty(i))-1.d0/xtv0220 ) ) xes = 1.0d0 elseif (ib.eq.14) then xlower = exp( dble(ee*nu12_1000) * @ ( 1.d0/dble(ty(i))-1.d0/xtv1000 ) ) xnu_trans = dble( nu(1,4)-nu12_1000 ) xtv_trans = xnu_trans / dble(nu(1,4)/tvt1(i)- @ nu12_1000/xtv1000) xes = (1.d0-exp( dble(-ee*xnu_trans/xtv_trans) )) / @ (1.d0-exp( dble(-ee*xnu_trans/ty(i)) )) elseif (ib.eq.13) then xlower = exp( dble(ee*nu12_0200) * @ ( 1.d0/dble(ty(i))-1.d0/xtv0200 ) ) xnu_trans = dble(nu(1,4)-nu12_0200) xtv_trans = xnu_trans / dble(nu(1,4)/tvt1(i)- @ nu12_0200/xtv0200) xes = (1.d0-exp( dble(-ee*xnu_trans/xtv_trans) )) / @ (1.d0-exp( dble(-ee*xnu_trans/ty(i)) )) else xlower = exp( dble(ee*elow(isot,ib)) * @ ( 1.d0/dble(ty(i))-1.d0/tvt0(i) ) ) xes = 1.0d0 end if xqv = (1.d0-exp( dble(-ee*667.3801/tvtbs(i)) )) / @ (1.d0-exp( dble(-ee*667.3801/ty(i)) )) xfactor = xlower * xqv**2.d0 * xes elseif (isot.eq.2) then xalfa = 1.d0/2.d0* exp( dble(-ee*(nu22_1000-nu(2,2))/ @ ty(i)) ) xbeta = 1.d0/2.d0* exp( dble(-ee*(nu22_0200-nu(2,2))/ @ ty(i)) ) xtv0200 = dble( - ee * nu22_0200 ) / @ ( log( xbeta/(1.d0+xalfa+xbeta) ) - dble(ee*nu(2,2))/ @ tvt0(i) ) xtv1000 = dble( - ee * nu22_1000 ) / @ ( log( xalfa/(1.d0+xalfa+xbeta) ) - dble(ee*nu(2,2))/ @ tvt0(i) ) if (ib.eq.14) then xlower = exp( dble(ee*nu22_1000) * @ ( 1.d0/dble(ty(i))-1.d0/xtv1000 ) ) xnu_trans = dble(nu(2,4)-nu22_1000) xtv_trans = xnu_trans / dble(nu(2,4)/tvt1(i)-nu22_1000/ @ xtv1000) xes = (1.d0-exp( dble(-ee*xnu_trans/xtv_trans) )) / @ (1.d0-exp( dble(-ee*xnu_trans/ty(i)) )) elseif (ib.eq.13) then xlower = exp( dble(ee*nu22_0200) * @ ( 1.d0/dble(ty(i))-1.d0/xtv0200 ) ) xnu_trans = dble( nu(2,4)-nu22_0200 ) xtv_trans = xnu_trans / dble(nu(2,4)/tvt1(i)-nu22_0200/ @ xtv0200) xes = (1.d0-exp( dble(-ee*xnu_trans/xtv_trans) )) / @ (1.d0-exp( dble(-ee*xnu_trans/ty(i)) )) else xlower = exp( dble(ee*elow(isot,ib)) * @ ( 1.d0/dble(ty(i))-1.d0/tvt0(i) ) ) xes = 1.0d0 end if xqv = (1.d0-exp( dble(-ee*662.3734/tvtbs(i)) )) / @ (1.d0-exp( dble(-ee*662.3734/ty(i)) )) xfactor = xlower * xqv**2.d0 * xes elseif (isot.eq.3) then xalfa = 1.d0/2.d0* exp( dble(-ee*(nu32_1000-nu(3,2))/ @ ty(i)) ) xbeta = 1.d0/2.d0* exp( dble(-ee*(nu32_0200-nu(3,2))/ @ ty(i)) ) xtv0200 = dble( - ee * nu32_0200 ) / @ ( log( xbeta/(1.d0+xalfa+xbeta) ) - dble(ee*nu(3,2))/ @ tvt0(i) ) xtv1000 = dble( - ee * nu32_1000 ) / @ ( log( xalfa/(1.d0+xalfa+xbeta) ) - dble(ee*nu(3,2))/ @ tvt0(i) ) if (ib.eq.14) then xlower = exp( dble(ee*nu32_1000) * @ ( 1.d0/dble(ty(i))-1.d0/xtv1000 ) ) xnu_trans = dble( nu(3,4)-nu32_1000 ) xtv_trans = xnu_trans / dble(nu(3,4)/tvt1(i)-nu32_1000/ @ xtv1000) xes = (1.d0-exp( dble(-ee*xnu_trans/xtv_trans) )) / @ (1.d0-exp( dble(-ee*xnu_trans/ty(i)) )) elseif (ib.eq.13) then xlower = exp( dble(ee*nu32_0200) * @ ( 1.d0/dble(ty(i))-1.d0/xtv0200 ) ) xnu_trans = dble( nu(3,4)-nu32_0200 ) xtv_trans = xnu_trans / dble(nu(3,4)/tvt1(i)-nu32_0200/ @ xtv0200) xes = (1.d0-exp( dble(-ee*xnu_trans/xtv_trans) )) / @ (1.d0-exp( dble(-ee*xnu_trans/ty(i)) )) else xlower = exp( dble(ee*elow(isot,ib)) * @ ( 1.d0/dble(ty(i))-1.d0/tvt0(i) ) ) xes = 1.0d0 end if xqv = (1.d0-exp( dble(-ee*648.4784/tvtbs(i)) )) / @ (1.d0-exp( dble(-ee*648.4784/ty(i)) )) xfactor = xlower * xqv**2.d0 * xes elseif (isot.eq.4) then xalfa = 1.d0/2.d0* exp( dble(-ee*(nu42_1000-nu(4,2))/ @ ty(i)) ) xbeta = 1.d0/2.d0* exp( dble(-ee*(nu42_0200-nu(4,2))/ @ ty(i)) ) xtv0200 = dble( - ee * nu42_0200 ) / @ ( log( xbeta/(1.d0+xalfa+xbeta) ) - dble(ee*nu(4,2))/ @ tvt0(i) ) xtv1000 = dble( - ee * nu42_1000 ) / @ ( log( xalfa/(1.d0+xalfa+xbeta) ) - dble(ee*nu(4,2))/ @ tvt0(i) ) if (ib.eq.14) then xlower = exp( dble(ee*nu42_1000) * @ ( 1.d0/dble(ty(i))-1.d0/xtv1000 ) ) xnu_trans = dble( nu(4,4)-nu42_1000 ) xtv_trans = xnu_trans / dble(nu(4,4)/tvt1(i)-nu42_1000/ @ xtv1000) xes = (1.d0-exp( dble(-ee*xnu_trans/xtv_trans) )) / @ (1.d0-exp( dble(-ee*xnu_trans/ty(i)) )) elseif (ib.eq.13) then xlower = exp( dble(ee*nu42_0200) * $ ( 1.d0/dble(ty(i))-1.d0/xtv0200 ) ) xnu_trans = dble( nu(4,4)-nu42_0200 ) xtv_trans = xnu_trans / dble(nu(4,4)/tvt1(i)-nu42_0200/ @ xtv0200) xes = (1.d0-exp( dble(-ee*xnu_trans/xtv_trans) )) / @ (1.d0-exp( dble(-ee*xnu_trans/ty(i)) )) else xlower = exp( dble(ee*elow(isot,ib)) * @ ( 1.d0/dble(ty(i))-1.d0/tvt0(i) ) ) xes = 1.0d0 end if xqv = (1.d0-exp( dble(-ee*664.7289/tvtbs(i)) )) / @ (1.d0-exp( dble(-ee*664.7289/ty(i)) )) xfactor = xlower * xqv**2.d0 * xes elseif (isot.eq.5 .and. ib.eq.1) then xlower = 1.d0 xes = 1.0d0 xqv = (1.d0-exp( dble(-ee*nuco_10/tvtbs(i)) )) / @ (1.d0-exp( dble(-ee*nuco_10/ty(i)) )) xfactor = xlower * xqv * xes end if con(i) = con(i) * xfactor if (i.eq.nzy) coninf = coninf * xfactor end do ! if ( isot.eq.1 .and. icurt_pop.eq.30 ) then ! close (30) ! endif return end