c***********************************************************
      
      SUBROUTINE gravi_out_gidl30(ichoice)
      
c***********************************************************
      INCLUDE 'ssbox.inc'                                       




c     nearby
      if(ichoice.eq.1) then

         gidl_version=3.
         write(30) gidl_version
         write(30) mx_mc,my_mc
         
         write(30) TIME,omega,lkm,bvak
         write(30) xl,yl,encmin,enclim
         write(30) nxg,nyg
         write(30) gmet,ngrav,roo
         write(30) nzg,fzrange
         
c     version 3

         write(30) apu_zmax,apu_zmean,apu_zdisp
         
         write(30) (s(i),i=1,lkm)
         write(30) (r(i,1),i=1,lkm)
         write(30) (r(i,2),i=1,lkm)
         write(30) (r(i,3),i=1,lkm)
         write(30) (agrav(i,1)/omega2,i=1,lkm)
         write(30) (agrav(i,2)/omega2,i=1,lkm)
         write(30) (agrav(i,3)/omega2,i=1,lkm)
      endif






c     here AGRAV includes near+grid+sheet
c     AGRAV_sheet = separately that due to sheet

      if(ichoice.eq.2) then
         write(30) (agrav(i,1)/omega2,i=1,lkm)
         write(30) (agrav(i,2)/omega2,i=1,lkm)
         write(30) (agrav(i,3)/omega2,i=1,lkm)
         write(30) (agrav_sheet(i)/omega2,i=1,lkm)

         write(30) tgrid,igrav
         nnxg=nxg
         nnyg=nyg
         nnzg=nzg
         write(30) 
     +        (((fxgrid3(i,j,k)/omega2,k=1,nnzg),j=1,nnyg),i=1,nnxg)
         write(30) 
     +        (((fygrid3(i,j,k)/omega2,k=1,nnzg),j=1,nnyg),i=1,nnxg)
         write(30) 
     +        (((fzgrid3(i,j,k)/omega2,k=1,nnzg),j=1,nnyg),i=1,nnxg)
         write(30) cpugrav1,cpugrav2,cpugrav3,cpugrav4
         
      endif

      end

c***********************************************************
      
      SUBROUTINE gravi_out_gidl31
      
c***********************************************************
      INCLUDE 'ssbox.inc'                                       

c     for tests
      real xtest(160),ytest(160),ztest(160)
      real fxapu31(160),fyapu31(160),fzapu31(160)
      

c     grid gravity at interpolated locations


c     from common
      zmin=apu_zmin
      zmax=apu_zmax



      n128x=nxg
      n128y=nyg
      n128z=2*nzg
      n64x=n128x/2
      n64y=n128y/2
      n64z=n128z/2
      dxg=2.d0*xl/n128x
      dyg=2.d0*yl/n128y
 
      n30=90
      write(31) time,omega,tgrid,igrav
      write(31) nxg,nyg,n30
      write(31) xl,yl,enclim,encmin
      
      do 61 i=1,20
         xtest(i)=-xl+0.05*(i-0.5)*2.*xl   
         ytest(i)=0.
         ztest(i)=0.
 61   continue
      do 62 i=1,10
         xtest(20+i)=0.
         ytest(20+i)=-yl+0.10*(i-0.5)*2.*yl   
         ztest(20+i)=0.
 62   continue
      do 621 i=1,60
         xtest(30+i)=0.
         ytest(30+i)=0.
         ztest(30+i)=zmin*2+(i-0.5)*(zmax-zmin)*4./60.
 621  continue

      do 63 i=1,n30
         
         x0p=xtest(i)
         y0p=ytest(i)-svak*x0p
         z0p=ztest(i)
         if(y0p.gt.yl)  y0p=y0p-2.*yl
         if(y0p.lt.-yl) y0p=y0p+2.*yl
         
         px=(x0p+xl)/2./xl*nxg+n128x-0.5
         py=(y0p+yl)/2./yl*nyg+n128y-0.5
         pz=(z0p-zmin)/(zmax-zmin)*nzg+n128z-0.5
         
         iix=int(px)-n128x
         iiy=int(py)-n128y
         iiz=int(pz)-n128z
         px=px-n128x
         py=py-n128y
         pz=pz-n128z
         
         dx=px-iix
         dy=py-iiy
         dz=pz-iiz
         tx=1.-dx
         ty=1.-dy
         tz=1.-dz
         
         ix1=iix+1
         ix2=iix+2
         iy1=iiy+1
         iy2=iiy+2
         iz1=iiz+1
         iz2=iiz+2
         
         if(ix1.gt.n128x) ix1=ix1-n128x
         if(ix2.gt.n128x) ix2=ix2-n128x
         if(iy1.gt.n128y) iy1=iy1-n128y
         if(iy2.gt.n128y) iy2=iy2-n128y
         if(ix1.lt.1) ix1=ix1+n128x
         if(ix2.lt.1) ix2=ix2+n128x
         if(iy1.lt.1) iy1=iy1+n128y
         if(iy2.lt.1) iy2=iy2+n128y

c     mass falling outside the grid
         if(iz2.gt.nzg) tz=0.
         if(iz2.lt.1) tz=0.
         if(iz1.gt.nzg) dz=0.
         if(iz1.lt.1) dz=0.
         if(iz1.lt.1) iz1=1
         if(iz2.lt.1) iz2=1
         if(iz1.gt.nzg) iz1=nzg
         if(iz2.gt.nzg) iz2=nzg

         fxapu=
     +        dx*dy*dz*fxgrid3(ix2,iy2,iz2)+
     +        dx*dy*tz*fxgrid3(ix2,iy2,iz1)+    
     +        dx*ty*dz*fxgrid3(ix2,iy1,iz2)+
     +        dx*ty*tz*fxgrid3(ix2,iy1,iz1)+     
     +        tx*dy*dz*fxgrid3(ix1,iy2,iz2)+
     +        tx*dy*tz*fxgrid3(ix1,iy2,iz1)+    
     +        tx*ty*dz*fxgrid3(ix1,iy1,iz2)+
     +        tx*ty*tz*fxgrid3(ix1,iy1,iz1)
         
         fyapu=
     +        dx*dy*dz*fygrid3(ix2,iy2,iz2)+
     +        dx*dy*tz*fygrid3(ix2,iy2,iz1)+    
     +        dx*ty*dz*fygrid3(ix2,iy1,iz2)+
     +        dx*ty*tz*fygrid3(ix2,iy1,iz1)+     
     +        tx*dy*dz*fygrid3(ix1,iy2,iz2)+
     +        tx*dy*tz*fygrid3(ix1,iy2,iz1)+    
     +        tx*ty*dz*fygrid3(ix1,iy1,iz2)+
     +        tx*ty*tz*fygrid3(ix1,iy1,iz1)

         fzapu=
     +        dx*dy*dz*fzgrid3(ix2,iy2,iz2)+
     +        dx*dy*tz*fzgrid3(ix2,iy2,iz1)+    
     +        dx*ty*dz*fzgrid3(ix2,iy1,iz2)+
     +        dx*ty*tz*fzgrid3(ix2,iy1,iz1)+     
     +        tx*dy*dz*fzgrid3(ix1,iy2,iz2)+
     +        tx*dy*tz*fzgrid3(ix1,iy2,iz1)+    
     +        tx*ty*dz*fzgrid3(ix1,iy1,iz2)+
     +        tx*ty*tz*fzgrid3(ix1,iy1,iz1)

c     _d2
c     if(abs(z0p).gt.zmax) fzapu=fzapu*abs(z0p)/zmax
c     _tracer
c     extrapolation: include forces beyond encmin 
c     encmin <  enclim

         if(abs(z0p).gt.zmax) then
            finf=2.*pii*gg*tmass/4./xl/yl
            fac=(1.+(encmin/z0p)**2)**(-.5)                            
            fzapu=-finf*fac*z0p/abs(z0p)
         endif
         if(abs(z0p).lt.zmax) then
            finf=2.*pii*gg*tmass/4./xl/yl
            fac=(1.+(enclim/z0p)**2)**(-.5)                            
            fzapu=fzapu-finf*fac*z0p/abs(z0p)
         endif

         fxapu31(i)=fxapu/omega2
         fyapu31(i)=fyapu/omega2
         fzapu31(i)=fzapu/omega2

 63   continue

      write(31) (xtest(i),i=1,n30)
      write(31) (ytest(i),i=1,n30)
      write(31) (ztest(i),i=1,n30)
      write(31) (fxapu31(i),i=1,n30)
      write(31) (fyapu31(i),i=1,n30)
      write(31) (fzapu31(i),i=1,n30)

      end




c***********************************************************
      
      SUBROUTINE gravi_out_terminal(ichoice)
      
c***********************************************************

      INCLUDE 'ssbox.inc'                                       


 1104 format(a10,f12.3,3e20.6)
 1105 format(a35,g14.3,2f12.3,i6)


c----------------------------------------------------------
c     first the contribution from near particles: fz_check1

      if(ichoice.eq.1) then

         fx_check1=0.
         fy_check1=0.
         fz_check1=0.
         
         axmin=1d8
         aymin=1d8
         azmin=1d8
         axmax=-1d8
         aymax=-1d8
         azmax=-1d8

         do 101 i=1,lkm
            fx_check1=fx_check1+xmass(i)*agrav(i,1)
            fy_check1=fy_check1+xmass(i)*agrav(i,2)
            fz_check1=fz_check1+xmass(i)*agrav(i,3)
            if(agrav(i,1).lt.axmin) axmin=agrav(i,1)
            if(agrav(i,2).lt.aymin) aymin=agrav(i,2)
            if(agrav(i,3).lt.azmin) azmin=agrav(i,3)
            if(agrav(i,1).gt.axmax) axmax=agrav(i,1)
            if(agrav(i,2).gt.aymax) aymax=agrav(i,2)
            if(agrav(i,3).gt.azmax) azmax=agrav(i,3)
 101     continue
         fx_check1=fx_check1/tmass/omega2
         fy_check1=fy_check1/tmass/omega2
         fz_check1=fz_check1/tmass/omega2
         
         axmin1=axmin/omega2
         aymin1=aymin/omega2
         azmin1=azmin/omega2
         axmax1=axmax/omega2
         aymax1=aymax/omega2
         azmax1=azmax/omega2

         write(6,*) '-------------------------------------------'
         write(6,1104) 'f_check1: ',orb,fx_check1,fy_check1,fz_check1
         write(6,1104) 'f_amin1:  ',orb,axmin1,aymin1,azmin1
         write(6,1104) 'f_amax1:  ',orb,axmax1,aymax1,azmax1
         write(6,1104) 'z-mean,disp:',orb,apu_zmean,apu_zdisp
      endif


c---------------------------------------------------------------
c     contribution from near+far particles: fz_check2

      if(ichoice.eq.2) then

            fx_check2=0.
            fy_check2=0.
            fz_check2=0.
            axmin=1d8
            aymin=1d8
            azmin=1d8
            axmax=-1d8
            aymax=-1d8
            azmax=-1d8
            n_zover=0
            do 102 i=1,lkm
               fx_check2=fx_check2+xmass(i)*agrav(i,1)
               fy_check2=fy_check2+xmass(i)*agrav(i,2)
               fz_check2=fz_check2+xmass(i)*agrav(i,3)
               if(agrav(i,1).lt.axmin) axmin=agrav(i,1)
               if(agrav(i,2).lt.aymin) aymin=agrav(i,2)
               if(agrav(i,3).lt.azmin) azmin=agrav(i,3)
               if(agrav(i,1).gt.axmax) axmax=agrav(i,1)
               if(agrav(i,2).gt.aymax) aymax=agrav(i,2)
               if(agrav(i,3).gt.azmax) azmax=agrav(i,3)
               if(abs(r(i,3)).gt.zmax) n_zover=n_zover

 102        continue
            fx_check2=fx_check2/tmass/omega2
            fy_check2=fy_check2/tmass/omega2
            fz_check2=fz_check2/tmass/omega2
            axmin2=axmin/omega2
            aymin2=aymin/omega2
            azmin2=azmin/omega2
            axmax2=axmax/omega2
            aymax2=aymax/omega2
            azmax2=azmax/omega2

            write(6,1104) 'f_check2:',orb,fx_check2,fy_check2,fz_check2
            write(6,1105) 'f_check2: zmean, zdisp, zmax, nov',
     +           apu_zmean,apu_zdisp,zmax,n_zover
            write(6,1104) 'f_amin2: ',orb,axmin2,aymin2,azmin2
            write(6,1104) 'f_amax2: ',orb,axmax2,aymax2,azmax2

         endif

c----------------------------------------------------------
c     contribution from near+far particles +
c     optional infinite sheet  (fzgrav3)


         if(ichoice.eq.3) then
            fz_check3=0.
            azmin=1d8
            azmax=-1d8
            apuapu=0.
            do 103 i=1,lkm
               fz_check3=fz_check3+xmass(i)*agrav(i,3)
               if(agrav(i,3).lt.azmin) azmin=agrav(i,3)
               if(agrav(i,3).gt.azmax) azmax=agrav(i,3)
               
 103        continue
            fz_check3=fz_check3/tmass/omega2
            azmin3=azmin/omega2
            azmax3=azmax/omega2
            write(6,1104) 'f_check3:',orb,fz_check3
            write(6,1104) 'f_amin3: ',orb,azmin3
            write(6,1104) 'f_amax3: ',orb,azmax3
         endif

c----------------------------------------------------------
c     add checks: sum of vertical gravity should be zero
c     optional gcorr

         if(ichoice.eq.4) then

            fz_check4=0.
            azmin=1d8
            azmax=-1d8
            do 104 i=1,lkm
               fz_check4=fz_check4+xmass(i)*agrav(i,3)
               if(agrav(i,3).lt.azmin) azmin=agrav(i,3)
               if(agrav(i,3).gt.azmax) azmax=agrav(i,3)
 104        continue
            fz_check4=fz_check4/tmass/omega2
            azmin4=azmin/omega2
            azmax4=azmax/omega2
            write(6,1104) 'f_check4:',orb,fz_check4
            write(6,1104) 'f_amin4: ',orb,azmin4
            write(6,1104) 'f_amax4: ',orb,azmax4
            write(6,*) '-------------------------------------------'

         endif

         end
