      subroutine density(isw)
      implicit real*8 (a-h,o-z)
      implicit integer*4 (i-n)
      save

      include "global.f"

      real*8 time,dose,weight,conc,adummy,tlag
      real*8 data
      real*8 Cl,V,expz1,expz2
      real*8 Q,QwCl,QwV
      real*8 A,Alag,E
      real*8 AwQ,AwQlag
      real*8 P,PwCl,PwV,PwS,Plag,PwClag,PwVlag,PwSlag
      real*8 O,S,f,fwO,fwS

      integer*4 i,j,isw
      integer*4 id0,id0lag,apgar,mdv,evid,iptr
      integer*4 count,point

      character*15 z1out, z2out

      integer*4 mrows,mcols,debug

      parameter (mrows=5,mcols=800,debug=0)

      dimension count(mcols),point(mcols)
      dimension data(mrows,mcols)


*     If isw=1, read and setup local data; compute nobs.

      if (isw.eq.1) then

        if (iunit9.eq.1) then
          j=0
          do i=1,12
            if (outfil(i:i).eq.' ') then
              j=i
              go to 5
            end if  
          end do
5         continue          
          if (j.eq.0) then
            z1out=outfil//'.z1'
            z2out=outfil//'.z2'
          else
            z1out=outfil(1:j-1)//'.z1'  
            z2out=outfil(1:j-1)//'.z2'  
          end if
          open (unit=21,file=z1out,status='unknown',form='formatted')
          open (unit=22,file=z2out,status='unknown',form='formatted')
        end if

        open (unit=14,file='/home/arg/r/mix/pheno/pheno.dat',
     &    status='old',form='formatted')

        nobs=0
        iptr=1
        count(iptr)=0
        point(iptr)=1
        id0lag=1

        do i=1,mcols
          read(14,14001,end=10)
     &      id0,time,dose,weight,apgar,conc,mdv,evid

          data(1,i)=time
          data(2,i)=dose
          data(3,i)=weight
          data(4,i)=conc
          if (apgar.lt.5.d0) then
            data(5,i)=1.d0
          else
            data(5,i)=0.d0
          end if

          if (conc.gt.0.d0) nobs=nobs+1

          if (id0lag.eq.id0) then
            count(iptr)=count(iptr)+1
          else
            iptr=iptr+1
            count(iptr)=1
            id0lag=id0
            point(iptr)=i
          end if
        end do
10      continue

        if (iptr.ne.ncases) then
          write(3,3001)
          stop
        end if

        close(14)
        return
      end if
 

*     If isw=2, compute the density.

      if (isw.eq.2) then
        weight=data(3,point(id))
        adummy=data(5,point(id))

        x(1)=weight
        x(2)=adummy

        expz1=dexp(z(1))
        expz2=dexp(z(2))

        Cl =  weight*tau(1)*expz1
        V  =  weight*tau(2)*(1.d0+adummy*tau(3))*expz2

        Q    = Cl/V
        QwCl = 1.d0/V
        QwV  = -Cl/(V *V)

        Plag   = 1.d0
        PwClag = 0.d0
        PwVlag = 0.d0
        PwSlag = 0.d0

        Alag   = 0.d0
        AwQlag = 0.d0
        tlag   = 0.d0
 
        do i=point(id),point(id)-1+count(id)
          time=data(1,i)
          dose=data(2,i)
          conc=data(4,i)
          E = dexp(-Q*(time-tlag))
          if (dose .gt. 0.d0) then
            A = Alag*E + dose
            AwQ = -Alag*(time-tlag)*E + AwQlag*E
          else
            A = Alag*E
            AwQ = -Alag*(time-tlag)*E + AwQlag*E
            O = A/V
            S = tau(4)
            call normalh(conc,O,S,2.d0,f,fwO,fwS)
            P   = Plag*f
            PwS  = Plag*fwS + PwSlag*f
            PwCl = Plag*fwO*(AwQ/V)*QwCl + PwClag*f
            PwV  = Plag*fwO*((AwQ/V)*QwV-A/(V*V)) + PwVlag*f
            Plag   = P
            PwSlag = PwS
            PwClag = PwCl
            PwVlag = PwV
          end if
          tlag=time
          Alag=A
          AwQlag=AwQ
        end do

        d = P
        dwtau(1)=PwCl*weight*expz1
        dwtau(2)=PwV*weight*(1.d0+adummy*tau(3))*expz2
        dwtau(3)=PwV*weight*tau(2)*adummy*expz2
        dwtau(4)=PwS
        dwz(1)=PwCl*Cl
        dwz(2)=PwV*V
        return
      end if

*     If isw=3, write individual effects to unit 9.

      if (isw.eq.3) then
        write(9,9001) id,Cl,V
        write(21,21000) z(1)
        write(22,22000) z(2)
      end if
        
*     If isw=4, print some debugging information.

      if (isw.eq.4) then

        if (iunit9.eq.1) then
          close(21)
          close(22)
        end if

        if ((iprint.eq.0).and.(debug.ne.0)) then
          do i=1,ncases
            write(3,3008) i,point(i),count(i)
            write(3,3009) (data(j,point(i)),j=1,mrows)
            write(3,3009) (data(j,point(i)-1+count(i)),j=1,mrows)
          end do
          write(3,3002) Cl,V
          write(3,3003) Q,QwCl,QwV
          write(3,3004) A,Alag,E
          write(3,3005) AwQ,AwQlag
          write(3,3006) P,PwCl,PwV,PwS,Plag,PwClag,PwVlag,PwSlag
          write(3,3007) O,S,f,fwO,fwS
          return
        end if
      end if

      return
3001  format(1x,'error, density, error reading data')
3002  format(1x,'Cl,V'/1x,2f20.8)
3003  format(1x,'R,RwCl,RwV'/1x,3f20.8)
3004  format(1x,'A,Alag,E'/1x,3f20.8)
3005  format(1x,'AwR,AwRlag'/1x,2f20.8)
3006  format(1x,'P,PwCl,PwV,PwS,Plag,PwClag,PwVlag,PwSlag'/1x,8f12.8)
3007  format(1x,'O,S,f,fwO,fwS'/1x,5f11.7)
3008  format(1x,'i,point,count',1x,3i5)
3009  format(1x,'data',1x,8f8.4)
3010  format(1x,'time, dose, A, conc, A/V',f7.1,4f10.6)
3011  format(1x,'time, f, fwO, fwS       ',f7.1,3f10.6)
3012  format(1x,'time, P, PwCl, PwV, PwS ',f7.1,4f10.6)
9001  format(i5,2d30.18)
14001 format(i7,3f7.0,i7,f7.0,i7,i2)
21000 format(1x,f28.20)
22000 format(1x,f28.20)
      end
