c*********************************************************************** c mzescape.f c*********************************************************************** c c program for calculating atmospheric escape functions, from c a calculation of transmittances and derivatives of these ones subroutine mzescape(ig,taustar,tauinf,tauii, ib,isot, iirw,iimu) c jul 2011 malv+fgg adapted to LMD-MGCM c nov 99 malv adapt mztf to compute taustar (pg.23b-ma c nov 98 malv allow for overlaping in the lorentz line c jan 98 malv version for mz1d. based on curtis/mztf.for c 17-jul-96 mlp&crs change the calculation of mr. c evitar: divide por cero. anhadiendo: ff c oct-92 malv correct s(t) dependence for all histogr bands c june-92 malv proper lower levels for laser bands c may-92 malv new temperature dependence for laser bands c @ 991 malv boxing for the averaged absorber amount and t c ? malv extension up to 200 km altitude in mars c 13-nov-86 mlp include the temperature weighted to match c the eqw in the strong doppler limit. c*********************************************************************** implicit none include 'nltedefs.h' include 'nlte_atm.h' include 'nlte_data.h' include 'nlte_curtis.h' include 'tcr_15um.h' include 'nlte_results.h' c arguments integer ig ! ADDED FOR TRACEBACK real*8 taustar(nl) ! o real*8 tauinf(nl) ! o real*8 tauii(nl) ! o integer ib ! i integer isot ! i integer iirw ! i integer iimu ! i c local variables and constants integer i, in, ir, im, k,j integer nmu parameter (nmu = 8) ! real*8 tauinf(nl) real*8 con(nzy), coninf real*8 c1, c2, ccc real*8 t1, t2 real*8 p1, p2 real*8 mr1, mr2 real*8 st1, st2 real*8 c1box(70), c2box(70) real*8 ff ! to avoid too small numbers real*8 tvtbs(nzy) real*8 st, beta, ts, eqwmu real*8 mu(nmu), amu(nmu) real*8 zld(nl), zyd(nzy) real*8 correc real deltanux ! width of vib-rot band (cm-1) character isotcode*2 real*8 maximum real*8 csL, psL, Desp, wsL ! for Strong Lorentz limit c formats 111 format(a1) 112 format(a2) 101 format(i1) 202 format(i2) 180 format(a80) 181 format(a80) c*********************************************************************** c some needed values ! rl=sqrt(log(2.d0)) ! pi2 = 3.14159265358989d0 beta = 1.8d0 ! imrco = 0.9865 c esto es para que las subroutines de mztfsub calculen we c de la forma apropiada para mztf, no para fot icls=icls_mztf c codigos para filenames ! if (isot .eq. 1) isotcode = '26' ! if (isot .eq. 2) isotcode = '28' ! if (isot .eq. 3) isotcode = '36' ! if (isot .eq. 4) isotcode = '27' ! if (isot .eq. 5) isotcode = '62' ! if(ib.eq.1.or.ib.eq.2.or.ib.eq.3.or.ib.eq.4.or.ib.eq.5 ! @ .or.ib.eq.6.or.ib.eq.7.or.ib.eq.8.or.ib.eq.9) then !! encode(2,101,ibcode1)ib ! write ( ibcode1, 101) ib ! else !! encode(2,202,ibcode2)ib ! write (ibcode2, 202) ib ! endif ! write (*,'( 30h calculating curtis matrix : ,2x, ! @ 8h band = ,i2,2x, 11h isotope = ,i2)') ib, isot c integration in angle !!!!!!!!!!!!!!!!!!!! c------- diffusivity approx. if (iimu.eq.1) then ! write (*,*) ' diffusivity approx. beta = ',beta mu(1) = 1.0d0 amu(1)= 1.0d0 c-------data for 8 points integration elseif (iimu.eq.4) then write (*,*)' 4 points for the gauss-legendre angle quadrature.' mu(1)=(1.0d0+0.339981043584856)/2.0d0 mu(2)=(1.0d0-0.339981043584856)/2.0d0 mu(3)=(1.0d0+0.861136311594053)/2.0d0 mu(4)=(1.0d0-0.861136311594053)/2.0d0 amu(1)=0.652145154862546 amu(2)=amu(1) amu(3)=0.347854845137454 amu(4)=amu(3) beta=1.0d0 c-------data for 8 points integration elseif(iimu.eq.8) then write (*,*)' 8 points for the gauss-legendre angle quadrature.' mu(1)=(1.0d0+0.183434642495650)/2.0d0 mu(2)=(1.0d0-0.183434642495650)/2.0d0 mu(3)=(1.0d0+0.525532409916329)/2.0d0 mu(4)=(1.0d0-0.525532409916329)/2.0d0 mu(5)=(1.0d0+0.796666477413627)/2.0d0 mu(6)=(1.0d0-0.796666477413627)/2.0d0 mu(7)=(1.0d0+0.960289856497536)/2.0d0 mu(8)=(1.0d0-0.960289856497536)/2.0d0 amu(1)=0.362683783378362 amu(2)=amu(1) amu(3)=0.313706645877887 amu(4)=amu(3) amu(5)=0.222381034453374 amu(6)=amu(5) amu(7)=0.101228536290376 amu(8)=amu(7) beta=1.0d0 end if c!!!!!!!!!!!!!!!!!!!!!!! ccc ccc determine abundances included in the absorber amount ccc c first, set up the grid ready for interpolation. do i=1,nzy zyd(i) = dble(zy(i)) enddo do i=1,nl zld(i) = dble(zl(i)) enddo c vibr. temp of the bending mode : 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 ) c 2nd: correccion a la n10(i) (cantidad de absorbente en el lower state) c por similitud a la que se hace en cza.for do i=1,nzy if (isot.eq.5) then con(i) = dble( coy(i) * imrco ) else con(i) = dble( co2y(i) * imr(isot) ) correc = 2.d0 * dexp( dble(-ee*elow(isot,2))/tvtbs(i) ) con(i) = con(i) * ( 1.d0 - correc ) endif c----------------------------------------------------------------------- c mlp & cristina. 17 july 1996 c change the calculation of mr. it is used for calculating partial press c alpha = alpha(self,co2)*pp +alpha(n2)*(pt-pp) c for an isotope, if mr is obtained by co2*imr(iso)/nt we are considerin c collisions with other co2 isotopes (including the major one, 626) c as if they were with n2. assuming mr as co2/nt, we consider collisions c of type 628-626 as of 626-626 instead of as 626-n2. c mrx(i)=con(i)/ntx(i) ! old malv ! mrx(i)= dble(co2x(i)/ntx(i)) ! mlp & crs c jan 98: c esta modif de mlp implica anular el correc (deberia revisar esto) mr(i) = dble(co2y(i)/nty(i)) ! malv, jan 98 c----------------------------------------------------------------------- end do ! como beta y 1.d5 son comunes a todas las weighted absorber amounts, ! los simplificamos: ! coninf = beta * 1.d5 * dble( con(n) / log( con(n-1) / con(n) ) ) coninf = dble( con(nzy) / log( con(nzy-1) / con(nzy) ) ) ! write (*,*) ' coninf =', coninf ccc ccc temp dependence of the band strength and ccc nlte correction factor for the absorber amount ccc call mztf_correccion ( coninf, con, ib, isot, 0 ) ccc ccc reads histogrammed spectral data (strength for lte and vmr=1) ccc !hfile1 = dirspec//'hi'//dn ! Ya no distinguimos entre d/n !! hfile1 = dirspec//'hid' ! (see why in his.for) ! hfile1='hid' !! if (ib.eq.13 .or. ib.eq.14 ) hfile1 = dirspec//'his' ! if (ib.eq.13 .or. ib.eq.14 ) hfile1 = 'his' ! if(ib.eq.1.or.ib.eq.2.or.ib.eq.3.or.ib.eq.4.or.ib.eq.5 ! @ .or.ib.eq.6.or.ib.eq.7.or.ib.eq.8.or.ib.eq.9) then ! if (isot.eq.1) hisfile = hfile1//'26-'//ibcode1//'.dat' ! if (isot.eq.2) hisfile = hfile1//'28-'//ibcode1//'.dat' ! if (isot.eq.3) hisfile = hfile1//'36-'//ibcode1//'.dat' ! if (isot.eq.4) hisfile = hfile1//'27-'//ibcode1//'.dat' ! if (isot.eq.5) hisfile = hfile1//'62-'//ibcode1//'.dat' ! else ! if (isot.eq.1) hisfile = hfile1//'26-'//ibcode2//'.dat' ! if (isot.eq.2) hisfile = hfile1//'28-'//ibcode2//'.dat' ! if (isot.eq.3) hisfile = hfile1//'36-'//ibcode2//'.dat' ! if (isot.eq.4) hisfile = hfile1//'27-'//ibcode2//'.dat' ! if (isot.eq.5) hisfile = hfile1//'62-'//ibcode2//'.dat' ! endif !write (*,*) ' /MZESCAPE/ hisfile: ', hisfile ! the argument to rhist is to make this compatible with mztf_comp.f, ! which is a useful modification of mztf.f (to change strengths of bands ! call rhist (1.0) if(ib.eq.1) then if(isot.eq.1) then !Case 1 mm=mm_c1 nbox=nbox_c1 tmin=tmin_c1 tmax=tmax_c1 do i=1,nbox_max no(i)=no_c1(i) dist(i)=dist_c1(i) do j=1,nhist sk1(j,i)=sk1_c1(j,i) xls1(j,i)=xls1_c1(j,i) xln1(j,i)=xln1_c1(j,i) xld1(j,i)=xld1_c1(j,i) enddo enddo do j=1,nhist thist(j)=thist_c1(j) enddo else if(isot.eq.2) then !Case 2 mm=mm_c2 nbox=nbox_c2 tmin=tmin_c2 tmax=tmax_c2 do i=1,nbox_max no(i)=no_c2(i) dist(i)=dist_c2(i) do j=1,nhist sk1(j,i)=sk1_c2(j,i) xls1(j,i)=xls1_c2(j,i) xln1(j,i)=xln1_c2(j,i) xld1(j,i)=xld1_c2(j,i) enddo enddo do j=1,nhist thist(j)=thist_c2(j) enddo else if(isot.eq.3) then !Case 3 mm=mm_c3 nbox=nbox_c3 tmin=tmin_c3 tmax=tmax_c3 do i=1,nbox_max no(i)=no_c3(i) dist(i)=dist_c3(i) do j=1,nhist sk1(j,i)=sk1_c3(j,i) xls1(j,i)=xls1_c3(j,i) xln1(j,i)=xln1_c3(j,i) xld1(j,i)=xld1_c3(j,i) enddo enddo do j=1,nhist thist(j)=thist_c3(j) enddo else if(isot.eq.4) then !Case 4 mm=mm_c4 nbox=nbox_c4 tmin=tmin_c4 tmax=tmax_c4 do i=1,nbox_max no(i)=no_c4(i) dist(i)=dist_c4(i) do j=1,nhist sk1(j,i)=sk1_c4(j,i) xls1(j,i)=xls1_c4(j,i) xln1(j,i)=xln1_c4(j,i) xld1(j,i)=xld1_c4(j,i) enddo enddo do j=1,nhist thist(j)=thist_c4(j) enddo else write(*,*)'isot must be 2,3 or 4 for ib=1!!' write(*,*)'stop at mzescape/312' stop endif else if (ib.eq.2) then if(isot.eq.1) then !Case 5 mm=mm_c5 nbox=nbox_c5 tmin=tmin_c5 tmax=tmax_c5 do i=1,nbox_max no(i)=no_c5(i) dist(i)=dist_c5(i) do j=1,nhist sk1(j,i)=sk1_c5(j,i) xls1(j,i)=xls1_c5(j,i) xln1(j,i)=xln1_c5(j,i) xld1(j,i)=xld1_c5(j,i) enddo enddo do j=1,nhist thist(j)=thist_c5(j) enddo else write(*,*)'isot must be 1 for ib=2!!' write(*,*)'stop at mzescape/336' stop endif else if (ib.eq.3) then if(isot.eq.1) then !Case 6 mm=mm_c6 nbox=nbox_c6 tmin=tmin_c6 tmax=tmax_c6 do i=1,nbox_max no(i)=no_c6(i) dist(i)=dist_c6(i) do j=1,nhist sk1(j,i)=sk1_c6(j,i) xls1(j,i)=xls1_c6(j,i) xln1(j,i)=xln1_c6(j,i) xld1(j,i)=xld1_c6(j,i) enddo enddo do j=1,nhist thist(j)=thist_c6(j) enddo else write(*,*)'isot must be 1 for ib=3!!' write(*,*)'stop at mzescape/360' stop endif else if (ib.eq.4) then if(isot.eq.1) then !Case 7 mm=mm_c7 nbox=nbox_c7 tmin=tmin_c7 tmax=tmax_c7 do i=1,nbox_max no(i)=no_c7(i) dist(i)=dist_c7(i) do j=1,nhist sk1(j,i)=sk1_c7(j,i) xls1(j,i)=xls1_c7(j,i) xln1(j,i)=xln1_c7(j,i) xld1(j,i)=xld1_c7(j,i) enddo enddo do j=1,nhist thist(j)=thist_c7(j) enddo else write(*,*)'isot must be 1 for ib=4!!' write(*,*)'stop at mzescape/384' stop endif else write(*,*)'ib must be 1,2,3 or 4!!' write(*,*)'stop at mzescape/389' endif if (isot.ne.5) deltanux = deltanu(isot,ib) if (isot.eq.5) deltanux = deltanuco c****** c****** calculation of tauinf(nl) c****** call initial ff=1.0e10 do i=nl,1,-1 if(i.eq.nl)then call intz (zl(i),c2,p2,mr2,t2, con) do kr=1,nbox ta(kr)=t2 end do ! write (*,*) ' i, t2 =', i, t2 call interstrength (st2,t2,ka,ta) aa = p2 * coninf * mr2 * (st2 * ff) bb = p2 * coninf * st2 cc = coninf * st2 dd = t2 * coninf * st2 do kr=1,nbox ccbox(kr) = coninf * ka(kr) ddbox(kr) = t2 * ccbox(kr) ! c2box(kr) = c2 * ka(kr) * beta * dble(deltaz) * 1.d5 c2box(kr) = c2 * ka(kr) * dble(deltaz) end do ! c2 = c2 * st2 * beta * dble(deltaz) * 1.d5 c2 = c2 * st2 * dble(deltaz) else call intz (zl(i),c1,p1,mr1,t1, con) do kr=1,nbox ta(kr)=t1 end do ! write (*,*) ' i, t1 =', i, t1 call interstrength (st1,t1,ka,ta) do kr=1,nbox ! c1box(kr) = c1 * ka(kr) * beta * dble(deltaz) * 1.d5 c1box(kr) = c1 * ka(kr) * dble(deltaz) end do ! c1 = c1 * st1 * beta * dble(deltaz) * 1.d5 c1 = c1 * st1 * dble(deltaz) aa = aa + ( p1*mr1*(c1*ff) + p2*mr2*(c2*ff)) / 2.d0 bb = bb + ( p1*c1 + p2*c2 ) / 2.d0 cc = cc + ( c1 + c2 ) / 2.d0 ccc = ( c1 + c2 ) / 2.d0 dd = dd + ( t1*c1 + t2*c2 ) / 2.d0 do kr=1,nbox ccbox(kr) = ccbox(kr) + @ ( c1box(kr) + c2box(kr) )/2.d0 ddbox(kr) = ddbox(kr) + @ ( t1*c1box(kr)+t2*c2box(kr) )/2.d0 end do mr2 = mr1 c2=c1 do kr=1,nbox c2box(kr) = c1box(kr) end do t2=t1 p2=p1 end if pt = bb / cc pp = aa / (cc*ff) ! ta=dd/cc ! tdop = ta ts = dd/cc do kr=1,nbox ta(kr) = ddbox(kr) / ccbox(kr) end do ! write (*,*) ' i, ts =', i, ts call interstrength(st,ts,ka,ta) ! call intershape(alsa,alna,alda,tdop) call intershape(alsa,alna,alda,ta) * ua = cc/st c next loop calculates the eqw for an especified path ua,pp,pt,ta eqwmu = 0.0d0 do im = 1,iimu eqw=0.0d0 do kr=1,nbox ua(kr) = ccbox(kr) / ka(kr) * beta * 1.0d5 / mu(im) if(ua(kr).lt.0.)write(*,*)'mzescape/480',ua(kr), $ ccbox(kr),ka(kr),beta,mu(im),kr,im,i,nl call findw (ig,iirw, 0, csL,psL, Desp, wsL) if ( i_supersat .eq. 0 ) then eqw=eqw+no(kr)*w else eqw = w + (no(kr)-1) * ( asat_box*dist(kr) ) endif end do eqwmu = eqwmu + eqw * mu(im)*amu(im) end do ! tauinf(i) = exp( - eqwmu / dble(deltanux) ) tauinf(i) = 1.d0 - eqwmu / dble(deltanux) if (tauinf(i).lt.0.d0) tauinf(i) = 0.0d0 if (i.eq.nl) then taustar(i) = 0.0d0 else taustar(i) = dble(deltanux) * (tauinf(i+1)-tauinf(i)) ! ~ / ( beta * cc * 1.d5 ) ~ / ( beta * ccc * 1.d5 ) endif end do ! i continue c****** c****** calculation of tau(in,ir) for n<=r c****** do 1 in=1,nl-1 call initial call intz (zl(in), c1,p1,mr1,t1, con) do kr=1,nbox ta(kr) = t1 end do call interstrength (st1,t1,ka,ta) do kr=1,nbox c1box(kr) = c1 * ka(kr) * dble(deltaz) end do c1 = c1 * st1 * dble(deltaz) call intz (zl(in+1), c2,p2,mr2,t2, con) do kr=1,nbox ta(kr) = t2 end do call interstrength (st2,t2,ka,ta) do kr=1,nbox c2box(kr) = c2 * ka(kr) * dble(deltaz) end do c2 = c2 * st2 * dble(deltaz) aa = aa + ( p1*mr1*(c1*ff) + p2*mr2*(c2*ff)) / 2.d0 bb = bb + ( p1*c1 + p2*c2 ) / 2.d0 cc = cc + ( c1 + c2 ) / 2.d0 dd = dd + ( t1*c1 + t2*c2 ) / 2.d0 do kr=1,nbox ccbox(kr) = ccbox(kr) + (c1box(kr)+c2box(kr))/2.d0 ddbox(kr) = ddbox(kr) + (t1*c1box(kr)+t2*c2box(kr))/2.d0 end do mr1=mr2 t1=t2 c1=c2 p1=p2 do kr=1,nbox c1box(kr) = c2box(kr) end do pt = bb / cc pp = aa / (cc * ff) ts = dd/cc do kr=1,nbox ta(kr) = ddbox(kr) / ccbox(kr) end do call interstrength(st,ts,ka,ta) call intershape(alsa,alna,alda,ta) eqwmu = 0.0d0 do im = 1,iimu eqw=0.0d0 do kr=1,nbox ua(kr) = ccbox(kr) / ka(kr) * beta * 1.0d5 / mu(im) if(ua(kr).lt.0.)write(*,*)'mzescape/566',ua(kr), $ ccbox(kr),ka(kr),beta,mu(im),kr,im,i,nl call findw (ig,iirw, 0, csL,psL, Desp, wsL) if ( i_supersat .eq. 0 ) then eqw=eqw+no(kr)*w else eqw = w + (no(kr)-1) * ( asat_box*dist(kr) ) endif end do eqwmu = eqwmu + eqw * mu(im)*amu(im) end do tauii(in) = exp( - eqwmu / dble(deltanux) ) !write (*,*) 'i,tauii=',in,tauii(in) 1 continue tauii(nl) = 1.0d0 c end return end