c*********************************************************************** c File with all subroutines required by mztf c c jan 98 malv basado en mztfsub_solar c jul 2011 malv+fgg adapted to LMD-MGCM c c contiene: c initial c intershape c interstrength c intz c rhist c interh c we c simrul c fi c f c findw c voigtf c*********************************************************************** c **************************************************************** subroutine initial c ma & crs !evita troubles 16-july-96 c **************************************************************** implicit none include 'nltedefs.h' include 'nlte_curtis.h' c local variables integer i c *************** eqw = 0.0d00 aa = 0.0d00 bb = 0.0d00 cc = 0.0d00 dd = 0.0d00 do i=1,nbox ua(i) = 0.0d0 ccbox(i) = 0.0d0 ddbox(i) = 0.0d0 end do return end c ********************************************************************** subroutine intershape(alsx,alnx,adx,xtemp) c interpolates the line shape parameters at a temperature xtemp from c input histogram data. c ********************************************************************** implicit none include 'nltedefs.h' include 'nlte_curtis.h' c arguments real*8 alsx(nbox_max),alnx(nbox_max),adx(nbox_max),xtemp(nbox_max) c local variables integer i, k c *********** ! write (*,*) 'intershape xtemp =', xtemp do 1, k=1,nbox if (xtemp(k).gt.tmax) then write (*,*) ' WARNING ! Tpath,tmax= ',xtemp(k),tmax xtemp(k) = tmax endif if (xtemp(k).lt.tmin) then write (*,*) ' WARNING ! Tpath,tmin= ',xtemp(k),tmin xtemp(k) = tmin endif i = 1 do while (i.le.mm) i = i + 1 if (abs(xtemp(k)-thist(i)) .lt. 1.0d-4) then !evita troubles alsx(k)=xls1(i,k) !16-july-1996 alnx(k)=xln1(i,k) adx(k)=xld1(i,k) goto 1 elseif ( thist(i) .le. xtemp(k) ) then alsx(k) = (( xls1(i,k)*(thist(i-1)-xtemp(k)) + @ xls1(i-1,k)*(xtemp(k)-thist(i)) )) / (thist(i-1)-thist(i)) alnx(k) = (( xln1(i,k)*(thist(i-1)-xtemp(k)) + @ xln1(i-1,k)*(xtemp(k)-thist(i)) )) / (thist(i-1)-thist(i)) adx(k) = (( xld1(i,k)*(thist(i-1)-xtemp(k)) + @ xld1(i-1,k)*(xtemp(k)-thist(i)) )) / (thist(i-1)-thist(i)) goto 1 end if end do write (*,*) @ ' error in xtemp(k). it should be between tmin and tmax' 1 continue return end c ********************************************************************** subroutine interstrength (stx, ts, sx, xtemp) c interpolates the line strength at a temperature xtemp from c input histogram data. c ********************************************************************** implicit none include 'nltedefs.h' include 'nlte_curtis.h' c arguments real*8 stx ! output, total band strength real*8 ts ! input, temp for stx real*8 sx(nbox_max) ! output, strength for each box real*8 xtemp(nbox_max) ! input, temp for sx c local variables integer i, k c *********** do 1, k=1,nbox ! if(xtemp(k).lt.ts)then ! write(*,*)'***********************' ! write(*,*)'mztfsub_overlap/EEEEEEH!',xtemp(k),ts,k ! write(*,*)'***********************' ! endif if (xtemp(k).gt.tmax) xtemp(k) = tmax if (xtemp(k).lt.tmin) xtemp(k) = tmin i = 1 do while (i.le.mm-1) i = i + 1 ! write(*,*)'mztfsub_overlap/136',i,xtemp(k),thist(i) if ( abs(xtemp(k)-thist(i)) .lt. 1.0d-4 ) then sx(k) = sk1(i,k) ! write(*,*)'mztfsub_overlap/139',sx(k),k,i goto 1 elseif ( thist(i) .le. xtemp(k) ) then sx(k) = ( sk1(i,k)*(thist(i-1)-xtemp(k)) + sk1(i-1,k)* @ (xtemp(k)-thist(i)) ) / (thist(i-1)-thist(i)) ! write(*,*)'mztfsub_overlap/144',sx(k),k,i goto 1 end if end do write (*,*) ' error in xtemp(kr) =', xtemp(k), @ '. it should be between ' write (*,*) ' tmin =',tmin, ' and tmax =',tmax stop 1 continue stx = 0.d0 if (ts.gt.tmax) ts = dble( tmax ) if (ts.lt.tmin) ts = dble( tmin ) i = 1 do while (i.le.mm-1) i = i + 1 ! write(*,*)'mztfsub_overlap/160',i,ts,thist(i) if ( abs(ts-thist(i)) .lt. 1.0d-4 ) then do k=1,nbox stx = stx + no(k) * sk1(i,k) ! write(*,*)'mztfsub_overlap/164',stx end do return elseif ( thist(i) .le. ts ) then do k=1,nbox stx = stx + no(k) * (( sk1(i,k)*(thist(i-1)-ts) + @ sk1(i-1,k)*(ts-thist(i)) )) / (thist(i-1)-thist(i)) ! write(*,*)'mztfsub_overlap/171',stx end do ! stop return end if end do return end c ********************************************************************** subroutine intz(h,aco2,ap,amr,at, con) c return interp. concentration, pressure,mixing ratio and temperature c for a input height h c ********************************************************************** implicit none include 'nltedefs.h' include 'nlte_atm.h' include 'nlte_curtis.h' c arguments real h ! i real*8 con(nzy) ! i real*8 aco2, ap, at, amr ! o c local variables integer k c ************ if ( ( h.lt.zy(1) ).and.( h.le.-1.e-5 ) ) then write (*,*) ' zp= ',h,' zy(1)= ',zy(1) stop'from intz: error in interpolation, z < minimum height' elseif (h.gt.zy(nzy)) then write (*,*) ' zp= ',h,' zy(nzy)= ',zy(nzy) stop'from intz: error in interpolation, z > maximum height' end if if (h.eq.zy(nzy)) then ap = dble( py(nzy) ) aco2= con(nzy) at = dble( ty(nzy) ) amr = dble( mr(nzy) ) return end if do k=1,nzy-1 if( abs( h-zy(k) ).le.( 1.e-5 ) ) then ap = dble( py(k) ) aco2= con(k) at = dble( ty(k) ) amr = dble( mr(k) ) return elseif(h.gt.zy(k).and.h.lt.zy(k+1))then ap = dble( exp( log(py(k)) + log(py(k+1)/py(k)) * @ (h-zy(k)) / (zy(k+1)-zy(k)) ) ) aco2 = exp( log(con(k)) + log( con(k+1)/con(k) ) * @ (h-zy(k)) / (zy(k+1)-zy(k)) ) at = dble( ty(k)+(ty(k+1)-ty(k))*(h-zy(k))/ @ (zy(k+1)-zy(k)) ) amr = dble( mr(k)+(mr(k+1)-mr(k))*(h-zy(k))/ @ (zy(k+1)-zy(k)) ) return end if end do 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 'nltedefs.h' include 'nlte_curtis.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 interh( sx, alsx, alnx, adx, xtemp, xtdop ) c interpolates a histogram at temperature xtemp from input histogr c jan 98 malv version para mz1d basada en el inth de solar10: c mz5/curtis/mztfsub_solar.f c **************************************************************** implicit none include 'nltedefs.h' include 'nlte_curtis.h' c arguments real*8 sx(nbox_max) ! o real*8 xtemp ! i real*8 alsx(nbox_max) ! o real*8 alnx(nbox_max) ! o real*8 adx(nbox_max) ! o real*8 xtdop ! i c local variables integer i, j, k c ************ if (xtemp.gt.thist(1)) then ! write (*,*) ' xtemp-path, thist(1)max: ', ! @ xtemp,thist(1) xtemp = thist(1) elseif (xtemp.lt.thist(nhist)) then ! write (*,*) ' xtemp-path, thist(nhist)min: ', ! @ xtemp,thist(nhist) xtemp = thist(nhist) end if i=0 1 i=i+1 if (abs(xtemp-thist(i)) .lt. 1.0d-4) goto 2 if (thist(i).lt.xtemp) goto 2 goto 1 2 j=i ! write (*,*) 'InterH/ j, xtemp, th(1),th(nh)', ! @ j, xtemp, thist(1),thist(nhist) if (j.gt.1) then do k=1,nbox sx(k) = ( sk1(j,k) * (thist(j-1)-xtemp) + sk1(j-1,k) * @ (xtemp-thist(j)) ) / (thist(j-1)-thist(j)) enddo elseif (j.eq.1) then do k=1,nbox sx(k) = sk1(1,k) enddo endif if (xtdop.gt.thist(1)) then ! write (*,*) ' xtdop-path, thist(1)max: ', ! @ xtdop,thist(1) xtdop = thist(1) elseif (xtdop.lt.thist(nhist)) then ! write (*,*) ' xtdop-path, thist(nhist)min: ', ! @ xtdop,thist(nhist) xtdop = thist(nhist) end if i=0 4 i=i+1 if (abs(xtdop-thist(i)) .lt. 1.0d-4) goto 5 if (thist(i).lt.xtdop) goto 5 goto 4 5 j=i ! write (*,*) 'InterH/ j, xtdop', ! @ j, xtdop if (j.gt.1) then do k=1,nbox alsx(k) = ( xls1(j,k) * (thist(j-1)-xtdop) + xls1(j-1,k)* @ (xtdop-thist(j)) ) / (thist(j-1)-thist(j)) alnx(k) = ( xln1(j,k) * (thist(j-1)-xtdop) + xln1(j-1,k)* @ (xtdop-thist(j)) ) / (thist(j-1)-thist(j)) adx(k) = ( xld1(j,k) * (thist(j-1)-xtdop) + xld1(j-1,k)* @ (xtdop-thist(j)) ) / (thist(j-1)-thist(j)) enddo elseif (j.eq.1) then do k=1,nbox alsx(k) = xls1(1,k) alnx(k) = xln1(1,k) adx(k) = xld1(j,k) enddo endif c end return end c ********************************************************************** real*8 function we(ig,me,pe,pl, idummy, nt_local,p_local, Desp, wsL) c icls=5 -->para mztf c icls=1,2,3-->para fot, line shape (v=1,l=2,d=3) (only use if wr=2) c calculates an approximate equivalent width for an error estimate. c c ioverlap = 0 ....... no correction for overlaping c 1 ....... "lisat" first correction (see overlap_box. c 2 ....... " " " plus "supersaturation" c idummy=0 do nothing c 1 write out some values for diagnostics c 2 correct the Strong Lorentz behaviour for SZA>90 c 3 casos 1 & 2 c malv nov-98 add overlaping's corrections c ********************************************************************** implicit none include 'nltedefs.h' include 'nlte_curtis.h' include 'tcr_15um.h' c arguments integer ig ! ADDED FOR TRACEBACK real*8 me ! I. path's absorber amount real*8 pe ! I. path's presion total real*8 pl ! I. path's partial pressure of CO2 real*8 nt_local ! I. needed for strong limit of Lorentz profil real*8 p_local ! I. " " " integer idummy ! I. indica varias opciones real*8 wsL ! O. need this for strong Lorentz correction real*8 Desp ! I. need this for strong Lorentz correction c local variables integer i real*8 y,x,wl,wd real*8 cn(0:7),dn(0:7) real*8 pi, xx real*8 f_sat_box real*8 dv_sat_box, dv_corte_box real*8 area_core_box, area_wing_box real*8 wlgood , parentesis , xlor real*8 wsl_grad c data blocks data cn/9.99998291698d-1,-3.53508187098d-1,9.60267807976d-2, @ -2.04969011013d-2,3.43927368627d-3,-4.27593051557d-4, @ 3.42209457833d-5,-1.28380804108d-6/ data dn/1.99999898289,5.774919878d-1,-5.05367549898d-1, @ 8.21896973657d-1,-2.5222672453,6.1007027481, @ -8.51001627836,4.6535116765/ c *********** c equivalent width of atmospheric line. pi = acos(-1.d0) if ( idummy.gt.9 ) @ write (*,*) ' S, m, alsa, pp =', ka(kr), me, alsa(kr), pl y=ka(kr)*me ! x=y/(2.0*pi*(alsa(kr)*pl+alna(kr)*(pe-pl))) x=y/(2.0d0*pi* alsa(kr)*pl) !+alna(kr)*(pe-pl))) ! Strong limit of Lorentz profile: WL = 2 SQRT( S * m * alsa*pl ) ! Para anular esto, comentar las siguientes 5 lineas ! if ( x .gt. 1.e6 ) then ! wl = 2.0*sqrt( y * alsa(kr)*pl ) ! else ! wl=y/sqrt(1.0d0+pi*x/2.0d0) ! endif wl=y/sqrt(1.0d0+pi*x/2.0d0) if (wl .le. 0.d0) then write(*,*)'mztfsub_overlap/496',ig,y,ka(kr),me,kr stop'WE/Lorentz EQW zero or negative!/498'!,ig endif if ( idummy.gt.9 ) @ write (*,*) ' y, x =', y, x xlor = x if ( (idummy.eq.2 .or. idummy.eq.12) .and. xlor.gt.1e5 ) then ! en caso que estemos en el regimen ! Strong Lorentz y la presion local ! vaya disminuyendo, corregimos la EQW ! con un gradiente analitico (notebook) wsL = 2.0*sqrt( y * alsa(kr)*pl ) wsl_grad = - 2.d0 * ka(kr)*alsa(kr) * nt_local*p_local / wsL wlgood = w_strongLor_prev(kr) + wsl_grad * Desp if (idummy.eq.12) @ write (*,*) ' W(wrong), W_SL, W_SL prev, W_SL corrected=', @ wl, wsL, w_strongLor_prev(kr), wlgood wl = wlgood endif ! wsL = wl pero esto no lo hacemos todavia, porque necesitamos ! el valor que ahora mismo tiene wsL para corregir la ! expresion R&W below ! write (*,*) 'WE arguments me,pe,pl =', me,pe,pl ! write (*,*) 'WE/ wl,ka(kr),alsa(kr) =', ! @ wl, ka(kr),alsa(kr) !>>>>>>> 500 format (a,i3,3(2x,1pe15.8)) 600 format (a,2(2x,1pe16.9)) 700 format (a,3(1x,1pe16.9)) ! if (kr.eq.8 .or. kr.eq.13) then ! write (*,500) 'WE/kr,m,pt,pl=', kr, me, pe, pl ! write (*,700) ' /aln,als,d_x=', alna(kr),alsa(kr), ! @ 2.0*pi*( alsa(kr)*pl + alna(kr)*(pe-pl) ) ! write (*,600) ' /alsa*p_CO2, alna*p_n2 :', ! @ alsa(kr)*pl, alna(kr)*(pe-pl) ! write (*,600) ' a*p, S =', ! @ alsa(kr)*pl + alna(kr)*(pe-pl) , ka(kr) ! write (*,600) ' /S*m, x =', y, x ! write (*,600) ' /aprox, WL =', ! @ 2.*sqrt( y*( alsa(kr)*pl+alna(kr)*(pe-pl) ) ), WL ! endif ! ! corrections to lorentz eqw due to overlaping and super-saturation ! i_supersat = 0 if ( icls.eq.5 .and. ioverlap.gt.0 ) then ! for the moment, only consider overlaping for mztf.f, not fot.f ! definition of saturation in the lisat model ! asat_box = 0.99d0 f_sat_box = 2.d0 * x xx = f_sat_box / log( 1./(1-asat_box) ) if ( xx .lt. 1.0d0 ) then dv_sat_box = 0.0d0 asat_box = 1.0d0 - exp( - f_sat_box ) else dv_sat_box = alsa(kr) * sqrt( xx - 1.0d0 ) ! approximation: only use of alsa in mars and venus endif ! area of saturated line ! area_core_box = 2.d0 * dv_sat_box * asat_box area_wing_box = 0.5d0 * ( wl - area_core_box ) dv_corte_box = dv_sat_box + 2.d0*area_wing_box/asat_box ! super-saturation or simple overlaping? ! ! i_supersat = 0 xx = dv_sat_box - ( 0.5d0 * dist(kr) ) if ( xx .ge. 0.0 ! definition of supersaturation @ .and. dv_sat_box .gt. 0.0 ! definition of saturation @ .and. (dist(kr).gt.0.0) ) ! box contains more than 1 line @ ! and not too far apart @ then i_supersat = 1 else ! no super-saturation, then use "lisat + first correction", i.e., ! correct for line products ! wl = wl endif end if ! end of overlaping loop if (icls.eq.2) then we = wl return endif cc doppler limit: if ( idummy.gt.9 ) @ write (*,*) ' S*m, alf_dop =', y, alda(kr)*sqrt(pi) x = y / (alda(kr)*sqrt(pi)) if ( x.lt.1.e-10 ) then ! to avoid underflow wd = y else wd=alda(kr)*sqrt(4.0*pi*x**2*(1.0+log(1.0+x))/(4.0+pi*x**2)) endif if ( idummy.gt.9 ) @ write (*,*) ' wd =', wd cc doppler weak limit c wd = ka(kr) * me cc good doppler if(icls.eq.5) then !para mztf !write (*,*) 'para mztf, icls=',icls if (x.lt.5.) then wd = 0.d0 do i=0,7 wd = wd + cn(i) * x**i end do wd = alda(kr) * x * sqrt(pi) * wd elseif (x.gt.5.) then wd = 0.d0 do i=0,7 wd = wd + dn(i) / (log(x))**i end do wd = alda(kr) * sqrt(log(x)) * wd else stop ' x should not be less than zero' end if end if if ( i_supersat .eq. 0 ) then parentesis = wl**2+wd**2-(wd*wl/y)**2 ! changed +(wd*wl/y)**2 to -...14-3-84 if ( parentesis .lt. 0.0 ) then if ((idummy.eq.2 .or. idummy.eq.12) .and. xlor.gt.1e5) then parentesis = wl**2+wd**2-(wd*wsL/y)**2 ! este cambio puede ser necesario cuando se hace ! correccion Strong Lor, para evitar valores ! negativos del parentesis en sqrt( ) else stop ' WE/ Error en las EQW wl,wl,y ' endif endif we = sqrt( parentesis ) ! write (*,*) ' from we: xdop,alda,wd', sngl(x),alda(kr),sngl(wd) ! write (*,*) ' from we: we', we else we = wl ! if there is supersaturation we can ignore wd completely; ! mztf.f will compute the eqw of the whole box afterwards endif if (icls.eq.3) we = wd if ( idummy.gt.9 ) @ write (*,*) ' wl,wd,w =', wl,wd,we wsL = wl return end c ********************************************************************** real*8 function simrul(a,b,fsim,c,acc) c adaptively integrates fsim from a to b, within the criterion acc. c ********************************************************************** implicit none real*8 res,a,b,g0,g1,g2,g3,g4,d,a0,a1,a2,h,x,acc,c,fsim real*8 s1(70),s2(70),s3(70) real*8 c1, c2 integer*4 m,n,j res=0. c=0. m=0 n=0 j=30 g0=fsim(a) g2=fsim((a+b)/2.) g4=fsim(b) a0=(b-a)*(g0+4.0*g2+g4)/2.0 1 d=2.0**n h=(b-a)/(4.0*d) x=a+(4.0*m+1.0)*h g1=fsim(x) g3=fsim(x+2.0*h) a1=h*(g0+4.0*g1+g2) a2=h*(g2+4.0*g3+g4) if ( abs(a1+a2-a0).gt.(acc/d)) goto 2 res=res+(16.0*(a1+a2)-a0)/45.0 m=m+1 c=a+m*(b-a)/d 6 if (m.eq.(2*(m/2))) goto 4 if ((m.ne.1).or.(n.ne.0)) goto 5 8 simrul=res return 2 m=2*m n=n+1 if (n.gt.j) goto 3 a0=a1 s1(n)=a2 s2(n)=g3 s3(n)=g4 g4=g2 g2=g1 goto 1 3 c1=c-(b-a)/d c2=c+(b-a)/d write(2,7) c1,c,c2,fsim(c1),fsim(c),fsim(c2) write(*,7) c1,c,c2,fsim(c1),fsim(c),fsim(c2) 7 format(2x,17hsimrule fails at ,/,3e15.6,/,3e15.6) goto 8 5 a0=s1(n) g0=g4 g2=s2(n) g4=s3(n) goto 1 4 m=m/2 n=n-1 goto 6 end c ********************************************************************** subroutine findw(ig,iirw,idummy,c1,p1, Desp, wsL) c this routine sets up accuracy criteria and calls simrule between limit c that depend on the number of atmospheric and cell paths. it gives eqw. c Add correction for EQW in Strong Lorentz regime and SZA>90 c ********************************************************************** implicit none include 'nltedefs.h' include 'nlte_curtis.h' c arguments integer ig ! ADDED FOR TRACEBACK integer iirw integer idummy ! I. indica varias opciones real*8 c1 ! I. needed for strong limit of Lorentz profil real*8 p1 ! I. " " " real*8 wsL ! O. need this for strong Lorentz correction real*8 Desp ! I. need this for strong Lorentz correction c local variables real*8 ept,eps,xa real*8 acc, c real*8 we real*8 f, fi, simrul external f,fi c ********** *********** ********* if(icls.eq.5) then !para mztf ! if(ig.eq.1682)write(*,*)'mztfsub_overlap/768',ua(kr),iirw if (iirw.eq.2) then !iirw=icf=2 ==> we use the w&r formula w = we(ig,ua(kr),pt,pp, idummy, c1,p1, Desp, wsL ) return end if ept=we(ig,ua(kr),pt,pp, idummy,c1,p1, Desp, wsL) else !para fot if (iirw.eq.2) then ! icf=2 ==> we use the w&r formula w = we(ig,sl_ua,pt,pp, idummy,c1,p1, Desp, wsL) return end if ept=we(ig,sl_ua,pt,pp, idummy,c1,p1, Desp, wsL) end if c the next block is a modification to avoid nul we. c this situation appears for weak lines and low path temperature, but c there is not any loss of accuracy. first july 1986 if (ept.eq.0.) then ! for weak lines sometimes we=0 ept=1.0e-18 write (*,*) 'ept =',ept write (*,*) 'from we: we=0.0' return end if acc = 4.d0 acc = 10.d0**(-acc) eps = acc * ept !accuracy 10-4 atmospheric eqw. xa=0.5*ept/f(0.d0) !width of doppler shifted atmospheric line. w=2.0*(simrul(0.0d0,xa,f,c,eps)+simrul(0.1d0,1.0/xa,fi,c,eps)) !no shift. return end c ********************************************************************** double precision function fi(y) c returns the value of f(1/y) c ********************************************************************** implicit none real*8 f, y fi=f(1.0/y)/y**2 return end c ********************************************************************** double precision function f(nu) c calculates 1-exp(-k(nu)u) for all series paths or combinations thereof c ********************************************************************** implicit none include 'nltedefs.h' include 'nlte_curtis.h' double precision tra,xa,ya,za,yy,nu double precision voigtf tra=1.0d0 yy=1.0d0/alda(kr) xa=nu*yy ya= ( alsa(kr)*pp + alna(kr)*(pt-pp) ) * yy za=ka(kr)*yy if(icls.eq.5) then !para mztf ! write (*,*) 'icls=',icls tra=za*ua(kr)*voigtf(sngl(xa),sngl(ya)) else tra=za*sl_ua*voigtf(sngl(xa),sngl(ya)) end if if (tra.gt.50.0) then tra=1.0 !2.0e-22 overflow cut-off. else if (tra.gt.1.0e-4) then tra=1.0-exp(-tra) end if f=tra return end c ********************************************************************** double precision function voigtf(x1,y) c computes voigt function for any value of x1 and any +ve value of y. c where possible uses modified lorentz and modified doppler approximatio c otherwise uses a rearranged rybicki routine. c c(n) = exp(-(n/h)**2)/(pi*sqrt(pi)), with h = 2.5 . c accurate to better than 1 in 10000. c ********************************************************************** implicit none real x1, y real x, xx, xxyy, xh,xhxh, yh,yhyh, f1,f2, p, q, xn,xnxn, voig real*8 b,g0,g1,g2,g3,g4,d1,d2,d3,d4,c integer*4 n dimension c(10) complex xp,xpp,z data c(1)/0.15303405/ data c(2)/0.94694928e-1/ data c(3)/0.42549174e-1/ data c(4)/0.13882935e-1/ data c(5)/0.32892528e-2/ data c(6)/0.56589906e-3/ data c(7)/0.70697890e-4/ data c(8)/0.64135678e-5/ data c(9)/0.42249221e-6/ data c(10)/0.20209868e-7/ x=abs(x1) if (x.gt.7.2) goto 1 if ((y+x*0.3).gt.5.4) goto 1 if (y.gt.0.01) goto 3 if (x.lt.2.1) goto 2 goto 3 c here uses modified lorentz approx. 1 xx=x*x xxyy=xx+y*y b=xx/xxyy voigtf=y*(1.+(2.*b-0.5+(0.75-(9.-12.*b)*b)/xxyy)/ * xxyy)/(xxyy*3.141592654) return c here uses modified doppler approx. 2 xx=x*x voigtf=0.56418958*exp(-xx)*(1.-y*(1.-0.5*y)*(1.1289-xx*(1.1623+ * xx*(0.080812+xx*(0.13854-xx*(0.033605-0.0073972*xx)))))) return c here uses a rearranged rybicki routine. 3 xh=2.5*x xhxh=xh*xh yh=2.5*y yhyh=yh*yh f1=xhxh+yhyh f2=f1-0.5*yhyh if (y.lt.0.1) goto 20 p=-y*7.8539816 !7.8539816=2.5*pi q=x*7.8539816 xpp=cmplx(p,q) z=cexp(xpp) d1=xh*aimag(z) d2=-d1 d3=yh*(1.-real(z)) d4=-d3+2.*yh voig=0.17958712*(d1+d3)/f1 goto 30 20 p=x*7.8539816 q=y*7.8539816 xp=cmplx(p,q) z=ccos(xp) d1=xh*aimag(z) d2=-d1 d3=yh*(1.-real(z)) d4=-d3+2.*yh voig=0.56418958*exp(y*y-x*x)*cos(2.*x*y)+0.17958712*(d1+d3)/f1 30 xn=0. do 55 n=1,10,2 xn=xn+1. xnxn=xn*xn g1=xh-xn g2=g1*(xh+xn) g3=0.5*g2*g2 voig=voig+c(n)*(d2*(g2+yhyh)+d4*(f1+xnxn))/(g3+yhyh*(f2+xnxn)) xn=xn+1. xnxn=xn*xn g1=xh-xn g2=g1*(xh+xn) g3=0.5*g2*g2 voig=voig+c(n+1)*(d1*(g2+yhyh)+d3*(f1+xnxn))/ @ (g3+yhyh*(f2+xnxn)) 55 continue voigtf=voig return end