subroutine hydrost(col,temp,zbot,rplnt,GMm,pres,z,den
     &   ,tmplst,natm)
C       --------------------------------------------------------------
C                    Hydrostatic Equilibrium Routine
C         Column density is the independent variable, P,Z calculated
C       --------------------------------------------------------------
	implicit double precision (a-h,o-z)
        include "parameters.dat"
	dimension pres(nd1),z(nd1),temp(nd1),col(nd1,0:4),den(nd1,0:4)
     &  , tmplst(nd1),rd(nd1)

C       -------------------------------------------------
C         Calcualte pressures at grid points
C              Start from top and integrate down
C       -------------------------------------------------
	rlast = zbot+rplnt
cd      plast=col(natm,0)*(GMm/rlast**2)
        plast = pres(natm)
	tlast = temp(natm)
        ulast = one/rlast
	z(natm)=zbot

	do n = natm-1, 1, -1
c   dcol = col(n,0) - col(n-1,0)
	   tnext = temp(n)
           pnext = pres(n)
	   tfac = half*rkb*(tnext+tlast)/GMm
           unext = ulast + tfac*dlog(pnext/plast)
           rnext = one/unext

cd   phalf1 = plast
cd 10        continue
cd      phalf0 = phalf1
cd      fac = one - half*thalf*dcol/rlast/phalf0
cd      pnext = plast + (GMm/rlast**2)*dcol/fac**2
cd      phalf1 = half*(pnext+plast)
cd      dp = dabs(phalf1-phalf0)
cd   if(dp .gt. 1.d-3) goto 10
cd         rnext = rlast - dcol*thalf/phalf1

           z(n) = rnext - rplnt
           tlast = tnext
           plast = pnext
           rlast = rnext
           ulast = unext
        enddo

C       --------------------------------------------------------
C            Calculate new densities, keep same mixing ratio
C       --------------------------------------------------------
        do n = 1, natm
           rd(n) = tmplst(n)/temp(n) 
        end do

        do i = 0, 4 
           do n = 1, natm 
              den(n,i) = rd(n)*den(n,i)
           end do
        end do

c        do i = 0, 4
c           col(1,i) = 
c           do n = 2, natm
c              col(n,i) = col(n-1,i)+half*(den(n,i)+den(n-1,i))
c     &        *(z(n-1)-z(n))
c           end do
c        end do

        return
        end