subroutine density(zt,prs,tmp,eddy,fmix,cp,den,col,ntau)
c       -------------------------------------------------------------
c                           Fixed parameters
c       -------------------------------------------------------------
        parameter (rkb=1.3807d-16,one=1.d0,half=5.d-1)

c       -------------------------------------------------------------
c                         External variables
c       -------------------------------------------------------------
        real*8 zt(1),prs(1),tmp(1),eddy,fmix(1),cp(1),den(65,0:2)
     &   ,col(65,0:2)
        integer*4 ntau

c       -------------------------------------------------------------
c                         Internal variables
c       ------------------------------------------------------------- 
        real*8 b1,b2,dum1,dum2,ht0
        integer*4 jt,im

c       ------------------------------------------------------------
        do jt = 1, ntau
           den(jt,0)=prs(jt)/rkb/tmp(jt)
           cp(jt) = 2.5d0*rkb*den(jt,0)
           b1 = 2.3d17*tmp(jt)**0.728d0
           b2 = 7.0d17*tmp(jt)**0.5d0
           dum1 = (one+b1/den(jt,0)/eddy)**-7
           dum2 = (one+b2/den(jt,0)/eddy)**-12
           den(jt,1)=dum2*fmix(1)*den(jt,0)
           den(jt,2)=dmin1(1.d9,dum2*fmix(2)*den(jt,0))

        enddo

 706    format(1x,1p5e12.4)

        jt = 1
        col(jt,1) =0.d0 
        col(jt,2) =0.d0 

        do jt = 2, ntau
           do im = 1, 2
              col(jt,im) = col(jt-1,im)+half*(den(jt,im)+den(jt-1,im))
     $         *(zt(jt-1)-zt(jt))
           enddo
           write(unit=*,fmt=706) zt(jt), prs(jt), col(jt,0), col(jt,1)
     &      , col(jt,2)

        enddo

        return
        end