c*********************************************************************** subroutine NLTE_leedat c reads planetary and molecular parameters c jan 98 malv first version c jul 2011 malv+fgg adapted to LMD-MGCM c*********************************************************************** implicit none include 'datafile.h' include 'nlte_paramdef.h' include 'nlte_commons.h' c local variables integer i,j, k,lun1, lun2 integer :: ib = 0 integer isot character isotcode*2 character ibcode1*1 c formats 132 format (i2) 101 format(i1) c*********************************************************************** lun1 = 1 lun2 = 2 do k=1,nisot ! encode (2,132,isotcode) indexisot(k) write (isotcode,132) indexisot(k) open (lun1, @ file=trim(datafile)//'/NLTEDAT/enelow'//isotcode//'.dat', @ status='old') open (lun2, @ file=trim(datafile)//'/NLTEDAT/deltanu'//isotcode//'.dat', @ status='old') read (lun1,*) read (lun2,*) read (lun1,*) (elow(k,i), i=1,nb) read (lun2,*) (deltanu(k,i), i=1,nb) close (lun1) close (lun2) end do c Call to rhist hfile1='hid' ! if (ib.eq.13 .or. ib.eq.14 ) hfile1 = dirspec//'his' if (ib.eq.13 .or. ib.eq.14 ) hfile1 = 'his' ib=1 do isot=1,4 c Cases 1,2,3,4 if (isot.eq.1) hisfile = hfile1//'26-1.dat' if (isot.eq.2) hisfile = hfile1//'28-1.dat' if (isot.eq.3) hisfile = hfile1//'36-1.dat' if (isot.eq.4) hisfile = hfile1//'27-1.dat' call rhist(1.0) if (isot.eq.1) then !Case 1 mm_c1=mm nbox_c1=nbox tmin_c1=tmin tmax_c1=tmax do i=1,nbox_max no_c1(i)=no(i) dist_c1(i)=dist(i) do j=1,nhist sk1_c1(j,i)=sk1(j,i) xls1_c1(j,i)=xls1(j,i) xln1_c1(j,i)=xln1(j,i) xld1_c1(j,i)=xld1(j,i) enddo enddo do j=1,nhist thist_c1(j)=thist(j) enddo else if(isot.eq.2) then !Case 2 mm_c2=mm nbox_c2=nbox tmin_c2=tmin tmax_c2=tmax do i=1,nbox_max no_c2(i)=no(i) dist_c2(i)=dist(i) do j=1,nhist sk1_c2(j,i)=sk1(j,i) xls1_c2(j,i)=xls1(j,i) xln1_c2(j,i)=xln1(j,i) xld1_c2(j,i)=xld1(j,i) enddo enddo do j=1,nhist thist_c2(j)=thist(j) enddo else if (isot.eq.3) then ! Case 3 mm_c3=mm nbox_c3=nbox tmin_c3=tmin tmax_c3=tmax do i=1,nbox_max no_c3(i)=no(i) dist_c3(i)=dist(i) do j=1,nhist sk1_c3(j,i)=sk1(j,i) xls1_c3(j,i)=xls1(j,i) xln1_c3(j,i)=xln1(j,i) xld1_c3(j,i)=xld1(j,i) enddo enddo do j=1,nhist thist_c3(j)=thist(j) enddo else if (isot.eq.4) then ! Case 4 mm_c4=mm nbox_c4=nbox tmin_c4=tmin tmax_c4=tmax do i=1,nbox_max no_c4(i)=no(i) dist_c4(i)=dist(i) do j=1,nhist sk1_c4(j,i)=sk1(j,i) xls1_c4(j,i)=xls1(j,i) xln1_c4(j,i)=xln1(j,i) xld1_c4(j,i)=xld1(j,i) enddo enddo do j=1,nhist thist_c4(j)=thist(j) enddo endif enddo do ib=2,4 isot=1 write (ibcode1,101) ib hisfile = hfile1//'26-'//ibcode1//'.dat' call rhist (1.0) if (ib.eq.2) then !Case 5 mm_c5=mm nbox_c5=nbox tmin_c5=tmin tmax_c5=tmax do i=1,nbox_max no_c5(i)=no(i) dist_c5(i)=dist(i) do j=1,nhist sk1_c5(j,i)=sk1(j,i) xls1_c5(j,i)=xls1(j,i) xln1_c5(j,i)=xln1(j,i) xld1_c5(j,i)=xld1(j,i) enddo enddo do j=1,nhist thist_c5(j)=thist(j) enddo else if (ib.eq.3) then !Case 6 mm_c6=mm nbox_c6=nbox tmin_c6=tmin tmax_c6=tmax do i=1,nbox_max no_c6(i)=no(i) dist_c6(i)=dist(i) do j=1,nhist sk1_c6(j,i)=sk1(j,i) xls1_c6(j,i)=xls1(j,i) xln1_c6(j,i)=xln1(j,i) xld1_c6(j,i)=xld1(j,i) enddo enddo do j=1,nhist thist_c6(j)=thist(j) enddo else if (ib.eq.4) then !Case 7 mm_c7=mm nbox_c7=nbox tmin_c7=tmin tmax_c7=tmax do i=1,nbox_max no_c7(i)=no(i) dist_c7(i)=dist(i) do j=1,nhist sk1_c7(j,i)=sk1(j,i) xls1_c7(j,i)=xls1(j,i) xln1_c7(j,i)=xln1(j,i) xld1_c7(j,i)=xld1(j,i) enddo enddo do j=1,nhist thist_c7(j)=thist(j) enddo endif enddo c end return end c ******************************************************************* subroutine rhist (factor_comp) c reads histogram data arrays created by ~/spectral/his.for c malv nov-98 add average distance between lines for overlapp c ******************************************************************* implicit none include 'nlte_paramdef.h' include 'nlte_commons.h' include 'datafile.h' c arguments real factor_comp c local variables integer j, r real*8 sk1_aux, xls1_aux, xln1_aux, xld1_aux,weight,nu0 character tonto*50 c formats ! 100 format(80a1) ! Esto es si fuese byte dummy(80) 100 format(a80) ! Esto es si fuese character dummy*80 150 format(a50) ! Esto es si fuese character dummy*80 c *************** open(unit=3, $ file=trim(datafile)//'/NLTEDAT/' $ //hisfile(1:len_trim(hisfile)),status='old') !read(3,100) dummy read(3,150) tonto read(3,*) weight read(3,*) mm read(3,*) nu0 read(3,*) nbox !read(3,'(a)') dumm read(3,'(a)') tonto if ( nbox .gt. nbox_max ) then write (*,*) ' nbox too large in input file ', hisfile stop ' Check maximum number nbox_max in mz1d.par ' endif do 1 j=1,mm ! for each temperature read(3,*) thist(j) do r=1,nbox ! for each box read(3,*) no(r), sk1(j,r), xls1(j,r),xln1(j,r),xld1(j,r), @ dist(r) c xld1(j,r)=xld1(j,r)*0.83255 !0.83255=sqrt(log(2)) enddo 1 continue tmax=thist(1) tmin=thist(mm) !close(unit=3,dispose='save') close(unit=3) do 2 j=1,mm do r=1,nbox sk1(j,r) = sk1(j,r) * factor_comp enddo 2 continue return end c*********************************************************************** subroutine leetvt c reads input vibr. temps. from external files or sets lte values c according to the driver table c jul 2011 malv+fgg adapted to LMD-MGCM c malv Jan 07 Add new vertical fine-grid for NLTE c jan 98 malv based on solar10sub c*********************************************************************** implicit none include 'nlte_paramdef.h' include 'nlte_commons.h' c local variables integer i real*8 zld(nl), zyd(nzy) real*8 xvt11(nzy), xvt21(nzy), xvt31(nzy), xvt41(nzy) c*********************************************************************** do i=1,nzy zyd(i) = dble(zy(i)) xvt11(i)= dble( ty(i) ) xvt21(i)= dble( ty(i) ) xvt31(i)= dble( ty(i) ) xvt41(i)= dble( ty(i) ) end do c interpolate to the nlte subgrid do i=1,nl zld(i) = dble( zl(i) ) enddo call interdp ( v626t1,zld,nl, xvt11,zyd,nzy, 1) call interdp ( v628t1,zld,nl, xvt21,zyd,nzy, 1) call interdp ( v636t1,zld,nl, xvt31,zyd,nzy, 1) call interdp ( v627t1,zld,nl, xvt41,zyd,nzy, 1) c end return end c*********************************************************************** subroutine getk (tt) c jan 98 malv version for mz1d. copied from solar10/getk.f c jul 2011 malv+fgg adapted to LMD-MGCM c*********************************************************************** implicit none include 'nlte_paramdef.h' include 'nlte_commons.h' c arguments real tt ! i. temperature !! local variables: real*8 k1x, k7x,k7xa,k7xb real*8 k3x, k3xaa,k3xac,k3xab,k3xbb,k3xba,k3xbc real*8 k3xco2,k3xn2,k3xco, k6x,k6x1,k6x2 real*8 k20x,k20xa,k20xb,k20xc real*8 k19xca,k19xcb,k19xcc real*8 k19xba,k19xbb,k19xbc real*8 k19xaa,k19xab,k19xac real*8 k21x,k21xa,k21xb,k21xc real*8 anu, factor , k21xa_626 real tt1,tt2, de integer i c*********************************************************************** co2i(001) + n2 ---> co2i + n2(1) k1 (not considered in the model c k1(i) i = 1 --- 626 ! 2 --- 636 ! 3 --- 628 ! 4 --- 627 ! k1x = 5.d-13 * sqrt(300.d0/tt) ! do i=1,nisot ! k1(i) = k1x * rf1 ! k1p(i) = k1(i) * exp( -ee/tt *1.d0* (-nun2+nu(i,4)) ) ! end do co2(001) + co2i ---> co2 + co2i(001) k2 (vv1 in table 4, paper i) c k2i i = x --- 628 ! y --- 636 ! z --- 627 k2a = 6.8d-12 * sqrt(tt) ! delta(e)< 42 cm-1 k2b = 3.6d-11 * sqrt(tt) * exp(-65.6557/26.3) ! > 42 cm-1 see table 3 k2a = k2a * rf2desac k2b = k2b * rf2desac k2x = 6.8d-12 * sqrt(tt) k2y = 3.6d-11 * sqrt(tt) * exp(-65.6557/26.3) k2z = 6.8d-12 * sqrt(tt) k2x = k2x * rf2iso k2y = k2y * rf2iso k2z = k2z * rf2iso k2xp = k2x * exp( dble( -ee/tt * (nu(1,4)-nu(2,4)) ) ) k2yp = k2y * exp( dble( -ee/tt * (nu(1,4)-nu(3,4)) ) ) k2zp = k2z * exp( dble( -ee/tt * (nu(1,4)-nu(4,4)) ) ) ! these are vt1 in table 3, paper i co2i(001) + m ---> co2i(0v0) + m k3 co2i(001) + o3p ---> co2i(0v0) + o3p k7 c k3vm(i) m = a --- co2 v = a --- 3 i = 1,2,3,4 ! b --- n2,o2 b --- 2 ! c --- co c k7v(i) v = a --- 3 i = 1,2,3,4 ! b --- 2 k7x = 2.0d-13*sqrt(tt/300.d0) k7x = k7x * rf7 if (iopt3.eq.0) then k3x = 2.2d-15 + 1.14d-10 * exp( dble(-76.75/tt**(1.d0/3.d0)) ) k3xaa = 3.0d-15 + 1.72d-10 * exp( dble(-76.75/tt**(1.d0/3.d0))) k3xac = 2.2d-15 + 9.66d-11 * exp( dble(-76.75/tt**(1.d0/3.d0))) k3x = k3x * rf3 k3xaa = k3xaa * rf3 k3xac = k3xac * rf3 tt1=220.0 !500.0 !220.0 tt2=250.0 !250.0 if(tt.le.tt1)then k3xab = k3x k3xbb = 0.d0 k7xa=k7x k7xb=0.d0 else if(tt.gt.tt2)then k3xab = 0.d0 k3xbb = k3x k7xa=0.d0 k7xb=k7x else k3xab = k3x/30.d0*(tt2-tt) k3xbb = k3x/30.d0*(tt-tt1) k7xa=k7x/30.d0*(tt2-tt) k7xb=k7x/30.d0*(tt-tt1) end if k3xba = k3xbb k3xbc = k3xbb elseif (iopt3.gt.0) then ! bauer et. al., 1987 if (tt.ge.190. .and. tt.le.250.) then factor = 0.9d0 + dble( (0.1-0.9)/(250.-190.) * (tt-190.) ) elseif (tt.lt.190.) then factor = 0.9d0 elseif (tt.gt.250.) then factor = 0.1d0 end if k3xn2 = 2.2d-15 + 1.14d-10 * exp(dble( -76.75/tt**(1.d0/3.d0))) k3xn2 = k3xn2 * rf3 k3xab = k3xn2 * factor k3xbb = k3xn2 * (1.-factor) k7xa = k7x * factor k7xb = k7x * (1.-factor) if (iopt3.eq.1) then if (tt.le.148.8) then k3xco2 = 13.8 - 0.1807 * (tt-140.) elseif (tt.ge.148.8 .and. tt.le.159.6) then k3xco2 = 12.21 - 0.1787 * (tt-148.8) elseif (tt.ge.159.6 .and. tt.le.171.0) then k3xco2 = 10.28 - 0.04035 * (tt-159.6) elseif (tt.ge.171.0 .and. tt.le.186.4) then k3xco2 = 9.82 - 0.027273 * (tt-171.0) elseif (tt.ge.186.4 .and. tt.le.244.1) then k3xco2 = 9.4 + 0.002253 * (tt-186.4) elseif (tt.ge.244.1 .and. tt.le.300) then k3xco2 = 9.53 + 0.02129 * (tt-244.1) elseif (tt.ge.300 .and. tt.le.336.1) then k3xco2 = 10.72 + 0.052632 * (tt-300) elseif (tt.ge.336.1 .and. tt.le.397.0) then k3xco2 = 12.62 + 0.0844 * (tt-336.1) elseif (tt.ge.397.0 .and. tt.le.523.4) then k3xco2 = 17.76 + 0.2615 * (tt-397.) end if k3xco2 = k3xco2 * 1.d-15 * rf3 k3xaa = 0.82d0 * k3xco2 k3xba = 0.18d0 * k3xco2 if (tt.le.163.) then k3xco = 10.58 - 0.093 * (tt-140) elseif (tt.ge.163. .and. tt.le.180.) then k3xco = 8.44 - 0.05353 * (tt-163.) elseif (tt.ge.180. .and. tt.le.196.) then k3xco = 7.53 - 0.034375 * (tt-180.) elseif (tt.ge.196. .and. tt.le.248.) then k3xco = 6.98 - 0.0108 * (tt-196.) elseif (tt.ge.248. .and. tt.le.301.) then k3xco = 6.42 + 0.01415 * (tt-248.) elseif (tt.ge.301. .and. tt.le.353.) then k3xco = 7.17 + 0.02865 * (tt-301.) end if k3xac = k3xco * 1.d-15 * rf3 k3xbc = 0.d0 elseif (iopt3.eq.2) then ! revision for the papers (feb 93) k3xco2 = 7.3d-14 * exp( -850.3/tt + 86523/tt**2.d0 ) k3xco2 = k3xco2 * rf3 k3xaa = 0.82d0 * k3xco2 k3xba = 0.18d0 * k3xco2 k3xco = 1.7d-14 * exp( -448.3/tt + 53636/tt**2.d0 ) k3xac = k3xco * rf3 k3xbc = 0.d0 end if end if do i=1,nisot k3aa(i) = k3xaa k3ba(i) = k3xba k3ab(i) = k3xab k3bb(i) = k3xbb k3ac(i) = k3xac k3bc(i) = k3xbc anu = nu(i,4)-nu(i,3) ! anu = nu(i,4)-nu(i,3)+70 k3aap(i) = k3aa(i) * exp( -ee/tt * anu )/6.d0 k3abp(i) = k3ab(i) * exp( -ee/tt * anu )/6.d0 k3acp(i) = k3ac(i) * exp( -ee/tt * anu )/6.d0 anu = nu(i,4)-nu(i,2) ! anu = nu(i,4)-nu(i,2)+40 k3bap(i) = k3ba(i) * exp( -ee/tt * anu )/4.d0 k3bbp(i) = k3bb(i) * exp( -ee/tt * anu )/4.d0 k3bcp(i) = k3bc(i) * exp( -ee/tt * anu )/4.d0 k7a(i) = k7xa k7b(i) = k7xb k7ap(i) = k7a(i) * exp(dble( -ee/tt*(nu(i,4)-nu(i,3)) ))/6.d0 k7bp(i) = k7b(i) * exp(dble( -ee/tt*(nu(i,4)-nu(i,2)) ))/4.d0 end do ! the next ones correspond to vv2 in table 4, paper i co2i(001) + co2 ---> co2i(020) + co2(010) k6 ! k6(i) i = 1,2,3,4 ! we need a new index for the inverse rates due to both fractions : c k6a(i) i=2,3,4 co2i(001) + co2 ---> co2i(020) + co2(010) c k6b(i) " co2i(001) + co2 ---> co2i(010) + co2(020) if (iopt6.eq.1) then if(tt.le.175.d0)then k6x=8.6d-15 elseif(tt.gt.175.0.and.tt.le.200.d0)then k6x=8.6d-15+9.d-16*(175.d0-tt)/25.d0 elseif(tt.gt.200.0.and.tt.le.225.d0)then k6x=7.7d-15+5.d-16*(200.d0-tt)/25.d0 elseif(tt.gt.225.0.and.tt.le.250.d0)then k6x=7.20d-15+6.d-16*(tt-225.d0)/25.d0 elseif(tt.gt.250.0.and.tt.le.275.d0)then k6x=7.80d-15+1.d-15*(tt-250.d0)/25.d0 elseif(tt.gt.275.0.and.tt.le.300.d0)then k6x=8.80d-15+1.3d-15*(tt-275.d0)/25.d0 elseif(tt.gt.300.0.and.tt.le.325.d0)then k6x=10.1d-15+1.54d-15*(tt-300.d0)/25.d0 elseif(tt.gt.325.0)then k6x=11.6d-15 end if elseif (iopt6.eq.2) then k6x = 3.6d-13 * exp( -1660/tt + 176948/tt**2.d0 ) if (tt.lt.175) k6x = 8.8d-15 end if k6x1 = k6x * rf6 * frac6 k6x2 = k6x * rf6 * (1.-frac6) k6 = k6x * rf6 k6p = k6 / 8.d0 * exp(dble( -ee/tt * (nu(1,4)-nu(1,2)-nu(1,1)) )) do i=2,nisot k6a(i) = k6x1 k6b(i) = k6x2 anu = nu(i,4)-nu(i,2)-nu(1,1) k6ap(i) = k6a(i) / 8.d0 * exp(dble( -ee*anu/tt )) anu = nu(i,4)-nu(i,1)-nu(1,2) k6bp(i) = k6b(i) / 8.d0 * exp(dble( -ee*anu/tt )) end do co2i(0v0) + co2i ---> co2i(0v-10) + co2i(010) c k5 esta reaccion es desdenable frente a co2 como colisionante. ! k5=3.0d-15+6.0d-17*(tt-210.d0) ! k5=k5*rf5 ! k5p=k5/2.d0*exp(-ee*125.77/tt) co2i(0v0) + m ---> co2i(0v-10) + m k19 (vt2,k5,k6 in table 3, paper co2i(0v0) + o3p ---> co2i(0v-10) + o3p k20 (vt2,k7 in table 3, paper c k19vm(i) m = a --- co2 v = a --- 3 i=1,2,3,4 ! b --- n2 b --- 2 ! c --- co c --- 1 c k20v(i) v = a --- 3 i = 1,2,3,4 ! b --- 2 ! c --- 1 ! ! k20x=1.9d-8*exp(-76.75/(tt**(1.d0/3.d0))) ! taylor,74 reajusted ! k20x=2.32d-9*exp(-76.75/(tt**(1./3.)))+1.0d-14*sqrt(tt) ! k&j, 83 ! k20x = 1.43d-12*(tt/300.d0) ! shved et al, 90 if (iopt20.eq.1) then ! first version of pap1 k20x=2.32d-9*exp(-76.75/(tt**(1./3.)))+3.5d-13*sqrt(tt) ! s&w, 91 k20xb = k20x / 2.d0 * rf20 k20xc = k20xb k20xa = 3.d0/2.d0 * k20xb elseif(iopt20.eq.2) then ! revision for the papers in feb 93 k20x=3.d-12 ! minimum value of lopez-puertas et al., 92 k20xc = k20x * rf20 k20xb = 2.d0 * k20xc k20xa = 3.d0/2.d0 * k20xb elseif(iopt20.eq.3) then ! values from boug & roble '91 k20x=1.d-12/sqrt(300.) * sqrt(tt) k20xc = k20x * rf20 k20xb = 2.d0 * k20xc k20xa = 3.d0/2.d0 * k20xb elseif(iopt20.eq.4) then ! values from boug & dick '88 case b k20x=7.d-13 k20xc = k20x * rf20 k20xb = 2.d0 * k20xc k20xa = 3.d0/2.d0 * k20xb elseif(iopt20.eq.5) then ! values from s.bougher (oct-98) k20x = 1.732d-13 * sqrt(tt) ! 1/sqrt(300) * sqrt(t) k20xc = k20x * rf20 k20xb = 2.d0 * k20xc k20xa = 3.d0/2.d0 * k20xb end if if (iopt19.eq.0) then k19xca = 4.64d-10 * exp(dble( - 74.75 / tt**(1.d0/3.d0) )) k19xcb = 6.69d-10 * exp(dble( - 84.07 / tt**(1.d0/3.d0) )) k19xcc = k19xcb if ( tt.le.250 ) then k19xba = 181.25d0 elseif ( tt.ge.310 ) then k19xba = 200.d0 + 0.9d0 * ( tt - 310.d0 ) else k19xba = 181.25d0 + 0.3125d0 * ( tt - 250.d0 ) end if k19xba = k19xba * 1.03558d-19 * tt ! cm-1 s-1 k19xbb = 1.24d-14 * ( tt / 273.3d0 )**2.d0 ! taine & lepoutre 1979 k19xbc = k19xbb k19xaa = 3.d0/2.d0 * k19xba k19xab = 3.d0/2.d0 * k19xbb k19xac = 3.d0/2.d0 * k19xbc k19xaa = k19xaa * rf19 k19xab = k19xab * rf19 k19xac = k19xac * rf19 k19xba = k19xba * rf19 k19xbb = k19xbb * rf19 k19xbc = k19xbc * rf19 k19xca = k19xca * rf19 k19xcb = k19xcb * rf19 k19xcc = k19xcc * rf19 elseif (iopt19.ge.1) then if (iopt19.eq.1) then ! lunt et. al., 1985 (thesis values) if (tt.le.175.) then k19xca = 4.d0 - 0.02d0 * (tt-140.d0) k19xcb = 0.494d0 + 0.0076 * (tt-140.d0) elseif (tt.ge.175. .and. tt.le.200.) then k19xca = 3.3d0 - 0.02d0 * (tt-175.) k19xcb = 0.76d0 + 0.0076d0 * (tt-175.d0) elseif (tt.ge.200. .and. tt.le.225.) then k19xca = 2.8d0 + 0.004d0 * (tt-200.d0) k19xcb = 0.95d0 + 0.014d0 * (tt-200.d0) elseif (tt.ge.225. .and. tt.le.250.) then k19xca = 2.9d0 + 0.024d0 * (tt-225.d0) k19xcb = 1.3d0 + 0.016d0 * (tt-225.d0) elseif (tt.ge.250. .and. tt.le.275.) then k19xca = 3.5d0 + 0.04d0 * (tt-250.d0) k19xcb = 1.7d0 + 0.032d0 * (tt-250.d0) elseif (tt.ge.275. .and. tt.le.295.) then k19xca = 4.5d0 + 0.055d0 * (tt-275.d0) k19xcb = 2.5d0 + 0.045d0 * (tt-275.d0) elseif (tt.ge.295. .and. tt.le.320.) then k19xca = 5.6d0 + 0.54d0 * (tt-295.d0) k19xcb = 3.4d0 + 0.045d0 * (tt-295.d0) end if k19xca = k19xca * 1.d-15 * rf19 k19xcb = k19xcb * 1.d-15 * rf19 k19xcc = k19xcb elseif (iopt19.eq.2) then ! revision for the papers, feb 1993 ! k19xca = 7.3d-14 * exp( -850.3d0/tt + 86523.d0/tt**2.d0 ) k19xca = 4.2d-12 * exp( -2988.d0/tt + 303930.d0/tt**2.d0 ) if (tt.le.175.) k19xca = 3.3d-15 k19xcb = 2.1d-12 * exp( -2659.d0/tt + 223052.d0/tt**2.d0 ) if (tt.le.175.) k19xcb = 7.6d-16 k19xca = k19xca * rf19 k19xcb = k19xcb * rf19 k19xcc = k19xcb elseif (iopt19.eq.3) then ! values from dick'72 for k19xc ! k19xcb is not modified if (tt.le.158.) then k19xca = 0.724d-15 elseif (tt.le.190.) then k19xca = 0.724d-15 + @ (1.1d-15-0.724d-15) * (tt-158.) / (190.-158.) elseif (tt.le.250.) then k19xca = 1.1d-15 + @ (3.45d-15-1.1d-15) * (tt-190.) / (250.-190.) elseif (tt.gt.250.) then k19xca = 3.45d-15 end if k19xcb = 2.1d-12 * exp( -2659.d0/tt + 223052.d0/tt**2.d0 ) if (tt.le.175.) k19xcb = 7.6d-16 k19xca = k19xca * rf19 k19xcb = k19xcb * rf19 k19xcc = k19xcb elseif (iopt19.eq.5) then k19xca = 5.2d-15 ! s.bougher, nov-98 k19xcb = 7.6d-16 ! nuestro, de feb-93 k19xcc = k19xcb k19xca = k19xca * rf19 k19xcb = k19xcb * rf19 end if factor = 2.5d0 k19xba = factor * k19xca k19xbb = factor * k19xcb k19xbc = factor * k19xcc factor = 3.d0/2.d0 k19xaa = factor * k19xba k19xab = factor * k19xbb k19xac = factor * k19xbc end if do i = 1, nisot k19aa(i) = k19xaa k19ba(i) = k19xba k19ca(i) = k19xca k19ab(i) = k19xab k19bb(i) = k19xbb k19cb(i) = k19xcb k19ac(i) = k19xac k19bc(i) = k19xbc k19cc(i) = k19xcc anu = nu(i,3)-nu(i,2) k19aap(i) = k19aa(i) * 6.d0/4.d0 * exp(dble( -ee*anu/tt)) k19abp(i) = k19ab(i) * 6.d0/4.d0 * exp(dble( -ee*anu/tt)) k19acp(i) = k19ac(i) * 6.d0/4.d0 * exp(dble( -ee*anu/tt)) anu = nu(i,2)-nu(i,1) k19bap(i) = k19ba(i) * 2.d0 * exp(dble( -ee*anu/tt)) k19bbp(i) = k19bb(i) * 2.d0 * exp(dble( -ee*anu/tt)) k19bcp(i) = k19bc(i) * 2.d0 * exp(dble( -ee*anu/tt)) anu = nu(i,1) k19cap(i) = k19ca(i) * 2.d0 * exp(dble( -ee*anu/tt)) k19cbp(i) = k19cb(i) * 2.d0 * exp(dble( -ee*anu/tt)) k19ccp(i) = k19cc(i) * 2.d0 * exp(dble( -ee*anu/tt)) k20a(i) = k20xa k20b(i) = k20xb k20c(i) = k20xc k20ap(i) = k20a(i)*6.d0/4.d0 * @ exp(dble( -ee/tt * (nu(i,3)-nu(i,2)) )) k20bp(i) = k20b(i)*4.d0/2.d0 * @ exp(dble( -ee/tt * (nu(i,2)-nu(i,1)) )) k20cp(i) = k20c(i)*2.d0/1.d0 * @ exp(dble( -ee/tt * nu(i,1) )) end do ! write(1,*) tt,k19cap(1),k19ac(1) ! the next ones correspond to vv3 in table 4 (paper i) co2i(0v0) + co2 ---> co2i(0v-10) + co2(010) k21 also see k33 c k21v(i) v = a --- 3 i = 1,2,3,4 ! b --- 2 ! c --- 1 ! we need a new index for the 030i rates due to both fractions : c k21a1 co2i(030) + co2 ---> co2i(020) + co2(010) c k21a2 co2i(030) + co2 ---> co2i(010) + co2(020) co2i(010) + co2j(000) ---> co2i(000) + co2j(010) kijk21c see pag.22-s ! k23k21c i=628,j=636 ! k24k21c i=628,j=627 ! k34k21c i=636,j=627 if (iopt21.eq.0) then k21x = 1.2d-11 k21xb = k21x k21xa = 3.d0/2.d0 * k21x k21xc = k21x ! esta ultima no se usa con 626 elseif (iopt21.eq.1) then k21x = 2.49d-11 ! orr & smith, 1987 k21xb = k21x k21xa = 3.d0/2.d0 * k21xb ! oscilador armonico k21xc = k21xb / 2.d0 ! novedad mia elseif (iopt21.eq.2) then k21x = 100.d0*k19xca ! dickinson'76 (icarus) k21xb = k21x k21xa = 3.d0/2.d0 * k21xb ! oscilador armonico k21xc = k21xb / 2.d0 ! novedad mia end if k21xa_626 = k21xa * rf21a !* 0.01d0 !* 10.d0 k21xa = k21xa * rf21a !* 0.01d0 k21xb = k21xb * rf21b k21xc = k21xc * rf21c k21a = k21xa_626 k21ap = k21a * 6.d0/8.d0 * @ exp( dble( -ee/tt * (nu(1,3)-nu(1,2)-nu(1,1)) ) ) do i = 2, nisot k21a1(i) = k21xa * frac21 k21a2(i) = k21xa * (1.d0-frac21) k21a1p(i) = k21a1(i) * 6.d0/8.d0 * @ exp(dble( -ee/tt* (nu(i,3)-nu(i,2)-nu(1,1)) )) k21a2p(i) = k21a2(i) * 6.d0/8.d0 * @ exp(dble( -ee/tt* (nu(i,3)-nu(i,1)-nu(1,2)) )) end do do i = 1, nisot k21b(i) = k21xb k21c(i) = k21xc k21bp(i) = k21b(i) * @ exp(dble( -ee/tt* (nu(i,2)-nu(i,1)-nu(1,1)) )) k21cp(i) = k21c(i) * @ exp(dble( -ee/tt * (nu(i,1)-nu(1,1)) )) end do k23k21c = k21xc k24k21c = k21xc k34k21c = k21xc k23k21cp = k23k21c*2.d0/2.d0 * @ exp(dble( -ee/tt* (nu(2,1)-nu(3,1)) )) k24k21cp = k24k21c*2.d0/2.d0 * @ exp(dble( -ee/tt* (nu(2,1)-nu(4,1)) )) k34k21cp = k34k21c*2.d0/2.d0 * @ exp(dble( -ee/tt* (nu(3,1)-nu(4,1)) )) ! these are also vv3 in table 4, paper i c k31 & k32 k31 = k21x * rf31 ! we're suposing thar the rate for the deactivation ! v-v from high combinational levels is the same k32 = k21x * rf32 ! that the one for : (020) --> (010) + (010) c o2(***) + co2i ---> co2(***) + co2i(***) k33 c k33a1 : co2(001) + co2i ---> co2(020) + co2i(010) (vv2, table 4, ! a2 : co2(001) + co2i ---> co2(010) + co2i(020) " ! b1 : co2(030) + co2i ---> co2(020) + co2i(010) (vv3, table 4, ! b2 : co2(030) + co2i ---> co2(010) + co2i(020) " ! c : co2(020) + co2i ---> co2(010) + co2i(010) " ! we have to add an index to the inverse rates, depending on the isotope k33c = k21x * rf33bc k33b1 = 3.d0/3.d0 * k33c * frac33 k33b2 = 3.d0/3.d0 * k33c * (1.d0-frac33) k33a1 = k6x * rf33a * frac33 k33a2 = k6x * rf33a * (1.d0-frac33) do i=2,nisot k33a1p(i)=k33a1* @ 1.d0/8.d0*exp(dble( -ee/tt* (nu(1,4)-nu(1,2)-nu(i,1)) )) k33a2p(i)=k33a2* @ 1.d0/8.d0*exp(dble( -ee/tt* (nu(1,4)-nu(1,1)-nu(i,2)) )) k33b1p(i)=k33b1* @ 6.d0/8.d0*exp(dble( -ee/tt* (nu(1,3)-nu(1,2)-nu(i,1)) )) k33b2p(i)=k33b2* @ 6.d0/8.d0*exp(dble( -ee/tt* (nu(1,4)-nu(1,1)-nu(i,2)) )) k33cp(i) =k33c * exp(dble( -ee/tt * (nu(1,2)-nu(1,1)-nu(i,1)) )) end do ! here they are the vt3 in table 3, paper i co2(2.7um) + m ---> co2(2.7um) + m k27 c k27a : n8 + m ---> n6 + m ! b : n7 + m ---> n6 + m ! c : n8 + m ---> n7 + m if (iopt27.eq.0) then k27a = 3.d-11 !between fermi levels k27b = 3.d-13 !between side levels k27c = 2.d0 * k27b !between side levels elseif (iopt27.ge.1) then ! orr & smith, 1987 k27a = 1.55d-12 k27c = 4.97d-12 k27b = k27c end if k27a = k27a * rf27f k27b = k27b * rf27s k27c = k27c * rf27s k27ap = k27a * exp(dble( -ee/tt * (nu(1,8)-nu(1,6)) )) k27bp = k27b * exp(dble( -ee/tt * (nu(1,7)-nu(1,6)) )) k27cp = k27c * exp(dble( -ee/tt * (nu(1,8)-nu(1,7)) )) ! the next two are not used in the model: c k28 : n* + n2 ---> n*low + n2(1) ! k28v v = a --- n8 ! b --- n7 ! c --- n6 ! k28a = 5.d-13 * sqrt(300.d0/tt) * rf28 ! = k1 ! k28b = k28a ! k28c = k28a ! k28ap = k28a * exp( -ee/tt * (nu(1,8)-1388.1847-nun2) ) ! k28bp = k28b * exp( -ee/tt * (nu(1,7)-1335.1317-nun2) ) ! k28cp = k28c * exp( -ee/tt * (nu(1,6)-1285.4087-nun2) ) c k29 : n* + co ---> n*low + co(1) ! k29v v = a --- n8 ! b --- n7 ! c --- n6 c k29a = ?????????? * rf29 c k29b = k29a c k29c = k29a c k29ap = k29a * exp( -ee/tt * (nu(1,8)-1388.1847-nuco) ) c k29bp = k29b * exp( -ee/tt * (nu(1,7)-1335.1317-nuco) ) c k29cp = k29c * exp( -ee/tt * (nu(1,6)-1285.4087-nuco) ) ! these are also vv1 processes in table 4, paper i c k26 : ! 1. deactivation of the 626 isotope: ! reaction: n* + co2i ---> n*low + co2i(001) ; i=1-4 ! nomenclature: k26v ; v=a,b,c,d for n8,n7,n6,n5 respectively ! inverse rates: k26vp(i) ; i=1-4 ! 2. deactivation of the minor isotopes: ! reaction: n*_i + co2j ---> n*_i_low + (001)_j ; i=2-4 ; j=1-4 ! nomenclature: k26vij ; v=a,c,d for n8,n6,n5 respectively ! inverse rates: k26vijp ! 3. notes: ! a * it is clear that : k26vij=/=k26vjip, ! b * at the moment we do not include inverse rates for the case 2., o ! the deactivation of the main isotope (pg. 32b, sn1). ! c * not 0221 level for minor isotopes is considered. ! d * only a value is known for these rates, so all of the deactivatio ! the same, but not the inverse rates. ! e * although all the direct deactivation constants have the same val ! it is useful to distinguish between them with the present nam k26a = 6.8d-12 * sqrt(tt) * rf26 ! = k2 k26b = k26a k26c = k26a if (iopt26.eq.0 .or. iopt26.eq.2) then k26d = k26a elseif( iopt26.eq.1) then k26d = 1.15d-10 * rf26 end if do i=1,4 k26ap(i) = k26a * @ exp(dble( -ee/tt * (nu(1,8)-nu12_1000-nu(i,4)) )) k26bp(i) = k26b * @ exp(dble( -ee/tt * (nu(1,7)- nu(1,2) -nu(i,4)) )) k26cp(i) = k26c * @ exp(dble( -ee/tt * (nu(1,6)-nu12_0200-nu(i,4)) )) k26dp(i) = k26d * @ exp(dble( -ee/tt * (nu(1,5)- nu(1,1) -nu(i,4)) )) end do k26a21 = k26a k26c21 = k26c k26d21 = k26d k26a22 = k26a k26c22 = k26c k26d22 = k26d k26a23 = k26a k26c23 = k26c k26d23 = k26d k26a24 = k26a k26c24 = k26c k26d24 = k26d k26a31 = k26a k26c31 = k26c k26d31 = k26d k26a32 = k26a k26c32 = k26c k26d32 = k26d k26a33 = k26a k26c33 = k26c k26d33 = k26d k26a34 = k26a k26c34 = k26c k26d34 = k26d k26a41 = k26a k26c41 = k26c k26d41 = k26d k26a42 = k26a k26c42 = k26c k26d42 = k26d k26a43 = k26a k26c43 = k26c k26d43 = k26d k26a44 = k26a k26c44 = k26c k26d44 = k26d !! some examples of inverse rates, although not used at the moment ! k26a32p = k26a32 * exp( -ee* (nu(3,8)-nu32_1000-nu(2,4)) / tt )*1./1. ! k26c32p = k26c32 * exp( -ee* (nu(3,6)-nu32_0200-nu(2,4)) / tt )*1./1. ! k26d32p = k26d32 * exp( -ee* (nu(3,5)- nu(3,1) -nu(2,4)) / tt )*2./2. ! ! k26a43 = k26a34 * exp( -ee* (nu(3,8)-nu32_1000-nu(4,4)) / tt )*1./1. ! k26c43 = k26c34 * exp( -ee* (nu(3,6)-nu32_0200-nu(4,4)) / tt )*1./1. ! k26d43 = k26d34 * exp( -ee* (nu(3,5)- nu(3,1) -nu(4,4)) / tt )*2./2. ! ! k26a24p = k26a24 * exp( -ee* (nu(2,8)-nu22_1000-nu(4,4)) / tt )*1./1. ! k26c24p = k26c24 * exp( -ee* (nu(2,6)-nu22_0200-nu(4,4)) / tt )*1./1. ! k26d24p = k26d24 * exp( -ee* (nu(2,5)- nu(2,1) -nu(4,4)) / tt )*2./2. ! this is taken as vv4 in table 4, paper i (in the inverse direction) c k41 : co(v) + co2 ---> co(v-1) + co2(001) + de ! k41_v v=1,2,3,4 ! ! de = -205.9 cm-1 when v=1 ! de = -232.9 cm-1 when v=2 ! de = -258.6 cm-1 when v=3 ! de = -285.0 cm-1 when v=4 k41p_taylor = 1.56d-11 * exp( -30.12/tt**0.333333 ) ! [ s-1 cm+3 ] k41p_shved = 7.5d7/sqrt(tt) ! [ s-1 atm-1 ] k41p_shved = k41p_shved * 1.38d-16/1013250. * tt ! [ s-1 cm+3 ] k41p_starr_hannock = 6.27d3 ! [ s-1 torr-1 ] if (iopt41.eq.1) then k41p_1 = k41p_starr_hannock * @ 760.*1.38d-16/1013250. * tt ! [ s-1 cm+3 ] elseif (iopt41.eq.2) then k41p_1 = 1.6d-12 * exp( -1169/tt + 77601/tt**2.d0 ) end if k41p_1 = k41p_1 * rf41 k41_1 = k41p_1 * exp(dble( -ee * 205.9/tt )) k41_2 = k41_1 k41p_2 = k41_2 * exp(dble( -ee * (-232.9)/tt )) k41_3 = k41_1 k41p_3 = k41_3 * exp(dble( -ee * (-258.6)/tt )) k41_4 = k41_1 k41p_4 = k41_4 * exp(dble( -ee * (-285.0)/tt )) !k41p_1 = k41p_1 * 1.d-6 !k41p_2 = k41p_2 * 1.d-6 c k41iso : 63(v) + co2 ---> 63(v-1) + co2(001) + de ! k41iso_v v=1,2,3 ! de = -253 cm-1 when v=1 ! de = -278 cm-1 when v=2 ! de = -303 cm-1 when v=3 k41iso_1 = k41_1 k41iso_1p = k41iso_1 * exp(dble( -ee * (-253.)/tt )) k41iso_2 = k41iso_1 k41iso_2p = k41iso_2 * exp(dble( -ee * (-278.)/tt )) k41iso_3 = k41iso_1 k41iso_3p = k41iso_3 * exp(dble( -ee * (-303.)/tt )) c k42 : co(v) + co ---> co(v-1) + co(1) + de=-26.481 si v=2 K42 ! -52.8940 si v=3 k42b ! -79.2402 si v=4 k42c ! tomado de stepanova & shved (ellos de powell, 1975), ver pg .. l5 ! solo para v=2 : k42 = 2.89d-10 * (1./sqrt(tt) + 67.4/tt**1.5) * @ exp(dble(24.78/tt)) k42 = k42 * rf42 k42b = k42 k42c = k42 k42p = k42 * exp(dble( -ee * (-26.481)/tt )) k42bp = k42b * exp(dble( -ee * (-52.894)/tt )) k42cp = k42c * exp(dble( -ee * (-79.24)/tt )) c k42iso : 63(v) + 63 ---> 63(v-1) + 63(1) + de=-25.31 si v=2 K42iso ! -50.57 si v=3 k42isob ! tomado de stepanova & shved (ellos de powell, 1975), ver pg .. l5 ! solo para v=2 : k42iso = k42 k42isop = k42iso * exp(dble( -ee * (-25.31)/tt )) k42isob = k42 k42isobp = k42isob * exp(dble( -ee * (-50.57)/tt )) c k43 : co(v) + o3p ---> co(v-1) + o3p + de=2143 ! tomado de lewittes et. al, 1978 para v=1 if (iopt43.eq.1) then tt1 = tt - 300. k43 = 2.85d-14 * exp( dble( 9.5e-3*tt1 + 1.11e-4*tt1**2. ) ) elseif (iopt43.eq.2) then k43 = 1.4d-5 * exp( -10957.d0 / tt + 1.486d6 / tt**2.d0 ) if ( tt.lt.265.0 ) k43 = 2.3d-14 end if k43 = k43 * rf43 k43p = k43 * exp( -ee * dble(2143.3 / tt) ) c k43iso : co63(v) + o3p ---> co63(v-1) + o3p + de=2096 ! Por similitud con el anterior k43iso = k43 k43isop = k43iso * exp( -ee * dble(2096. / tt) ) c k44 : co62(v) + co63 ---> co62(v-1) + co63(1) + de ! basado en Lopez-Valverde et al para el caso v=1, solo usamos este ! k44x x = a --- v=1 de= 147.33 ! b --- v=2 de= 20.7241 ! c --- v=3 de= -5.7 ! d --- v=4 de= -32.0361 k44a = 2.0d-12 * rf44 ! Solo vamos a usar este, no los b,c,d k44b = k44a k44c = k44a k44d = k44a de = 147.33 k44ap = k44a * exp(dble( -ee * de/tt )) de = 20.7241 k44bp = k44b * exp(dble( -ee * de/tt )) de = -5.7 k44cp = k44c * exp(dble( -ee * de/tt )) de = -32.0361 k44dp = k44d * exp(dble( -ee * de/tt )) co2(hcl) + co2 --> co2 + co2 + de(hcl) ! este rate tambien lo usamos para los high combination levels (para tra ! al lte. cualquier valor vale, supongo. es k_vthcl k_vthcl = 3.3d-15 ! similar al valor pequenho del vt2 k_vthcl = k_vthcl * rf_hcl return end