       program testmodels
       open(32,file='ringsneutc999',status='new')
       open(31,file='otrajc999',status='old')
c       open(32,file='testit4-rings',status='new')
c       open(31,file='trajdat3-d',status='old')
       tdum=0.
       istart=0
       iopt=2
c       iopt=1
c above is to read ascii file on first call
1      read(31,*,end=99)a,x,y,z
       r=sqrt(x**2+y**2+z**2)
       theta=acos(z/r)
       phi=atan2(y,x)
       theta=theta*45./atan(1.)
       phi=phi*45./atan(1.)
c insert model sub call here
       call satrings(r,theta,phi,waterden,istart,iopt)   
cc       if (r.gt.5) then
cc          write(32,*)x,y,z,tdum,0.0
cc       else
          write(32,*)x,y,z,a,waterden
cc       endif
c       write(32,*)x,y,z,tdum,waterden
       tdum=tdum+1.
       go to 1
99     stop
       end
      subroutine satrings(r,theta,phi,ringwat,istart,iopt)
c coordinates are Saturn polar coords
c densities in cm-3 of water for rings atmosphere or rings+moons
c iopt=1 is rings only
c      2=rings+all icy moons contributing
c
       dimension rarray(104,81)
       common/gridpts/xrd(104),xlatd(81)
       if(istart.gt.0)go to 10
       istart=1
       do 346 j=1,104
       do 346 k=1,81
       rarray(j,k)=-9.
346    continue
c set up grid values of r, lat
       do 347 jr=1,104
       xrd(jr)=(jr-1)*.25+(.25/2.)
347    continue
       do 348 ilat=1,81
       xlatd(ilat)=90.-ilat*(180./120.)+.75
348    continue
c read data for selected option
       imin=40
       imax=80
       if(iopt.ne.1)go to 245
       open(30,file='ringdata',status='old')
       go to 345
245    if(iopt.gt.2)return
       open(30,file='ringmoons',status='old')
345    continue
       read(30,300)((rarray(j,i),i=40,80),j=1,100)
300    format(6(6x,f7.4))   
       close(30)
10     continue
       uselat=90.-theta
       uselat=abs(uselat)
       user=r
c interpolate for densities
       call interprings(user,uselat,imin,imax,rarray,ringwat)
       ringwat=10.**ringwat
       if(ringwat.le.1.e-9)ringwat=0.
       return
       end
CC********************************************************
CC    Subroutine to define grid rectangle surrounding point
CC********************************************************
      subroutine interprings(user,uselat,imin,imax,rarray,quant)
      dimension xrd(104),xlatd(81),rarray(104,81)
      common/gridpts/xrd,xlatd
       quant=-9.
       if(user.gt.xrd(100).or.uselat.gt.xlatd(1))return
       if(user.lt.xrd(5).or.uselat.lt.0.)return
      do ix=1,104
       if(xrd(ix).le.user) jx=ix
      enddo
      do ix=1,61
       if(xlatd(ix).gt.uselat) jy=ix
      enddo
       call ringscalc(jx,jy,user,uselat,xrd,xlatd,rarray,quant)
      RETURN
      END
CC
CC*******************************************************************
CC  Subroutine to interpolate for density
CC*******************************************************************
      subroutine ringscalc(jx,jy,xpt,ypt,x,y,dw,quant)
CC
      dimension x(104),y(81),dw(104,81)
       quant=-9.
      kx=jx+1
      ky=jy+1
       if(kx.gt.104.or.ky.gt.61)return
      X9=(xpt-x(jx))/(x(kx)-x(jx))
      Y9=(ypt-y(jy))/(y(ky)-y(jy))
CC
      BLLL=dw(jx,jy) 
      BHLL=dw(kx,jy)
      BLHL=dw(jx,ky)
      BHHL=dw(kx,ky)
      VAL3 =BLHL+X9*(BHHL-BLHL)
      VAL4 =BLLL+X9*(BHLL-BLLL)
      VAL7 =BHLL+Y9*(BHHL-BHLL)
      VAL8 =BLLL+Y9*(BLHL-BLLL)
CC
      VB5=0.5*(VAL8 +X9*(VAL7 -VAL8 )+VAL4+Y9*(VAL3-VAL4))
CC
      quant=VB5
      RETURN
      END
c****************************************************************
