       program testmodels
       dimension isp(5),algn(5)
       character*64 fname
       open(28,file='testneutfiles',status='old')
c      open(29,file='mhd/flybysout',status='old')
       tdum=0.
       istart=0
       write(*,'("enter model, 1=Yung, 2=Fox/Yelle")')
       read(*,*)nmod
       write(*,'("3 species ids? 1-10 Yung, 1-40 Fox/Yelle")')
       read(*,*)isp(1),isp(2),isp(3)
c above is to read ascii file on first call
101    format(a)
2      continue
       open(32,file='oTAatm',status='old')
       read(28,101,end=99)fname
       open(33,file='neutyungTA10',status='new')
c       read(29,101,end=99)fname
1      read(32,*,end=98)x,y,z,tca
       r=sqrt(x**2+y**2+z**2)
       alt=(r-1.)*2575.
c insert model sub call here
c required coord system is TSO (Titan-centered, referenced to Sun)
c nmod=1 is Yung, nmod=2 is Yelle
       xno=1.e-33
       do 3 jspec=1,3
       ispecn=isp(jspec)
       call titneuts(r,x,nmod,ispecn,xnn,istart)   
       if(xnn.le.0.)xnn=xno
       algn(jspec)=alog10(xnn)
3      continue
       write(33,*)tca,(algn(i),i=1,3),alt
       go to 1
98     continue
       close(32)
       close(33)
c      go to 2
99     stop
       end
c***********************************************************
       subroutine titneuts(r,x,nmod,ispecn,xnn,istart)
c required coord system is TSO (Titan-centered, Sun referenced x axis)
c Titan neutral atmosphere models from Yung et al. and Yelle
c nmod=1 is Yung et al. model (ApJ Suppl.v55, p465, 1984)
c nmod=2 is Fox/Yelle model (personal communication, 1997) 
c also see GRL  24, p 2179, 1997)
c
c Species identifiers (ispecn) for Yung et al. model:
c 1=N2, 2=CH4, 3=H2, 4=HCN, 5=C2H2, 6=C2H4, 7=C2H6, 8=N
c 9=NH, 10=H
c altitude range is 725 to 9025 km
c
c Species identifiers for Fox/Yelle model:
c 1=N2,2=CH4,3=C2H2,4=C2H4,5=C2H6,6=C4H2,7=H,8=H2,9=C, 10=CH,
c 11=CH2,12=1CH2,13=CH3,14=C2H,15=C2H3,16=C2H5,17=C4H3,18=N,19=N2D,
c 20=NH,21=NH2,22=NH3,23=CN,24=HCN,25=H2CN,26=CNN,27=CHCN,28=C2N2,    
c 29=HC2N2,30=C3N,31=C2HCN,32=H2C3N,33=CH3CN,34=METAM,35=C4N2,
c 36=C2H3CN,37=C2H5CN,38=H2O,39=O,40=CO
c altitude range is 600-3000 km
c
c densities in cm-3
c
       character*64 fneuname
       common/modeln/dneut(10,831)             
       common/modeln2/yneut(41,161)
       if(istart.gt.0)go to 10
       istart=1
101    format(a)
       open(30,file='fneulist',status='old')
       read(30,101)fneuname
       open(31,file=fneuname,status='old')
       m=1
11     continue
       read(31,*,end=9)altd,(dneut(k,m),k=1,4)
       m=m+1
       go to 11
9      continue
       close(31)
       read(30,101)fneuname
       open(31,file=fneuname,status='old')
       m=1       
12     continue
       read(31,*,end=13)altd,(dneut(k,m),k=5,10)
       m=m+1
       go to 12
13     continue
       close(31)
       read(30,101)fneuname
       open(31,file=fneuname,status='old')
       do 16 i=1,41
       j=1
       do 17 m=1,20
       read(31,300)(yneut(i,k),k=j,j+7)
300    format(1p8e15.8)
       j=j+8
17     continue
       read(31,301)yneut(i,j)
301    format(1pe15.8)
16     continue
       close(31)
       close(30)
10     continue
c      compute altitude and solar zenith angle for interpolation
       alt=2575.*r-2575.
       xnn=0.
       if(nmod.eq.1.and.alt.le.725.)return
       if(nmod.eq.1.and.alt.ge.9025.)return
       if(nmod.eq.1)call fnmod(den,alt,ispecn)
       if(nmod.eq.2.and.alt.le.600.)return
       if(nmod.eq.2.and.alt.ge.3000.)return
       if(nmod.eq.2)call ynmod(den,alt,ispecn)
       if(nmod.lt.1.or.imod.gt.2)return
       xnn=den
       return
       end
C------------------------------------------------------------ 
      subroutine fnmod(den,alt,ISPECN)
c     this program calculates the neutral densities according to
c the model of Yung et al. 1984
      common/modeln/dneut(10,831)
      den=0.
      nalt=(alt-725.)/10.+1
       if(alt.lt.725.or.nalt.gt.830)go to 99
       stuff1=dneut(ISPECN,nalt+1)
       stuff2=dneut(ISPECN,nalt)
       den=stuff2+(alt-725.-(float(nalt)-1.)*10.)/10.*(stuff1-stuff2)
99     continue
       return
       end
c---------------------------------------------------------------------- 
      subroutine ynmod(den,alt,ISPECN)
c this subroutine calculates the Yelle neutral atmosphere
      common/modeln2/yneut(41,161)
      den=0.
      nalt=(alt-600.)/10.+1
      if(alt.lt.600.or.nalt.gt.160)go to 99
      stuff1=yneut(ISPECN,nalt+1)
      stuff2=yneut(ISPECN,nalt)
      den=stuff2+(alt-600.-(float(nalt)-1.)*10.)/10.*(stuff1-stuff2)
99    continue
      return
      end
