*...:....1....:....2....:....3....:....4....:....5....:....6....:....7....:....8
      subroutine density(isw)
      implicit real*8 (a-h,o-z)
      implicit integer*4 (i-n)
      save

      include "../code/global.f"
      include "local.f"

      real*8 delta,Fa,Fe,timeL
      real*8 O,OL,Owtau,OwtauL,Owz,OwzL
      real*8 Oa,OaL,Oawtau,OawtuL,Oawz,OawzL
      real*8 ka,kawtau,kawz
      real*8 ke,kewtau,kewz
      real*8 V,Vwtau,Vwz
      real*8 P,Pwtau,Pwz,PL,PwtauL,PwzL
      real*8 S,f,fwO,fwS
      real*8 OaVsav,Osave

      integer*4 i,j,isw
      integer*4 id0,id0lag,iptr
      integer*4 count,point,obs
      integer*4 idsave,ctsave

      character*12 fname

      dimension count(mcols),point(mcols),obs(mcols)
      dimension kawtau(mtau),kawz(mM)
      dimension kewtau(mtau),kewz(mM)
      dimension Vwtau(mtau),Vwz(mM)
      dimension Owtau(mtau),OwtauL(mtau),Owz(mM),OwzL(mM)
      dimension Oawtau(mtau),OawtuL(mtau),Oawz(mM),OawzL(mM)
      dimension Pwtau(mtau),Pwz(mM),PwtauL(mtau),PwzL(mM)
      dimension OaVsav(mcols),Osave(mcols)
      dimension idsave(mcols),ctsave(mcols)

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

      if (isw.eq.1) then

        if (iunit9.eq.1) then 
          open (unit=14,file='quinplot.dat',
     &      status='old',form='formatted')
        else
          open (unit=14,file='quin.dat',
     &      status='old',form='formatted')
        endif
        
        nobs=0
        iptr=1
        count(iptr)=0
        point(iptr)=1
        obs(iptr)=0
        id0lag=1

        do i=1,mcols

          ctsave(i) = 0

          read(14,14001,end=10)
     &      id0,
     &      time(i),conc(i),dose(i),
     &      form(i),observ(i),change(i),steady(i),
     &      dinter(i),age(i),height(i),weight(i),
     &      race(i),smoker(i),lush(i),
     &      heart(i),creatn(i),
     &      album(i),glyco(i)

          if (observ(i).eq.0) nobs=nobs+1

          if ((id0lag.eq.id0).or.(i.eq.1)) then
            count(iptr)=count(iptr)+1
            if (observ(i).eq.0) obs(iptr)=obs(iptr)+1
          else
            iptr=iptr+1
            count(iptr)=1
            point(iptr)=i
            if (observ(i).eq.0) then
              obs(iptr)=1
            else
              obs(iptr)=0
            end if
            id0lag=id0
          end if

          idsave(i) = iptr
          ctsave(i) = count(iptr)
          OaVsav(i) = 0.d0
          Osave(i)  = 0.d0

        end do

10      continue

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

        close(14)
        return
      end if


*     If isw=2, compute the density.

      if (isw.eq.2) then

        timeL=0.d0

        OaL=0.d0
        OL=0.d0
        do j=1,ltau
          Oawtau(j)=0.d0
          OawtuL(j)=0.d0
          Owtau(j)=0.d0
          OwtauL(j)=0.d0
        end do
        do j=1,M
          Oawz(j)=0.d0
          OawzL(j)=0.d0
          Owz(j)=0.d0
          OwzL(j)=0.d0
        end do

        PL = 1.d0
        do j=1,ltau
          Pwtau(j)=0.d0
          PwtauL(j)=0.d0
        end do
        do j=1,M
          Pwz(j)=0.d0
          PwzL(j)=0.d0
        end do

        do i=point(id),point(id)-1+count(id)


          index = i
          call g(ka,kawtau,kawz,ke,kewtau,kewz,V,Vwtau,Vwz)

*         Steady state or not?
          if (steady(i).ge.1) then

*           Steady state.

            Fa = 1.d0/(1.d0 - dexp(-ka*dinter(i)))
            Fe = 1.d0/(1.d0 - dexp(-ke*dinter(i)))
            Oa = (dose(i)/V)*Fa
            O  = (dose(i)/V)*(ka/(ka-ke))*(Fe-Fa)
            do j=1,ltau-1
              Oawtau(j) = 
     &          - (Oa/V)*Vwtau(j)
     &          - (dose(i)/V)*Fa**2*dinter(i)*dexp(-ka*dinter(i))
     &              *kawtau(j)
              Owtau(j) = 
     &          - (O/V)*Vwtau(j)
     &          + (dose(i)/V)/(ka-ke)*(Fe-Fa)*kawtau(j)
     &          - (dose(i)/V)*ka/(ka-ke)**2*(Fe-Fa)*kawtau(j)
     &          + (dose(i)/V)*ka/(ka-ke)*Fa**2
     &            *dinter(i)*dexp(-ka*dinter(i))*kawtau(j)
     &          + (dose(i)/V)*ka/(ka-ke)**2*(Fe-Fa)*kewtau(j)
     &          - (dose(i)/V)*ka/(ka-ke)*Fe**2
     &            *dinter(i)*dexp(-ke*dinter(i))*kewtau(j)
            end do
            do j=1,M
              Oawz(j) = 
     &          - (Oa/V)*Vwz(j)
     &          - (dose(i)/V)*Fa**2*dinter(i)*dexp(-ka*dinter(i))
     &              *kawz(j)
              Owz(j) = 
     &          - (O/V)*Vwz(j)
     &          + (dose(i)/V)/(ka-ke)*(Fe-Fa)*kawz(j)
     &          - (dose(i)/V)*ka/(ka-ke)**2*(Fe-Fa)*kawz(j)
     &          + (dose(i)/V)*ka/(ka-ke)*Fa**2
     &            *dinter(i)*dexp(-ka*dinter(i))*kawz(j)
     &          + (dose(i)/V)*ka/(ka-ke)**2*(Fe-Fa)*kewz(j)
     &          - (dose(i)/V)*ka/(ka-ke)*Fe**2
     &            *dinter(i)*dexp(-ke*dinter(i))*kewz(j)
            end do

          else

*           Non-steady state.

            delta = time(i)-timeL
            Fa = dexp(-ka*delta)
            Fe = dexp(-ke*delta)
            Oa = OaL*Fa+(dose(i)/V)
            O  = OL*Fe + OaL*(ka/(ka-ke))*(Fe-Fa)
            do j=1,ltau-1
              Oawtau(j) =
     &          + OawtuL(j)*Fa 
     &          - OaL*delta*Fa*kawtau(j) 
     &          - dose(i)*Vwtau(j)/(V**2)    
              Owtau(j) =
     &          + OwtauL(j)*Fe 
     &          - OL*delta*Fe*kewtau(j)
     &          + OawtuL(j)*ka*(Fe-Fa)/(ka-ke)
     &          + OaL*(Fe-Fa)*kawtau(j)/(ka-ke)
     &          + OaL*ka*(Fe-Fa)*(kewtau(j)-kawtau(j))/((ka-ke)**2)
     &          + OaL*ka*delta*(kawtau(j)*Fa-kewtau(j)*Fe)/(ka-ke)
            end do
            do j=1,M
              Oawz(j) =
     &            + OawzL(j)*Fa 
     &            - OaL*delta*Fa*kawz(j)
     &            - dose(i)*Vwz(j)/(V**2)    
              Owz(j) =
     &            + OwzL(j)*Fe 
     &            - OL*delta*Fe*kewz(j)
     &            + OawzL(j)*ka*(Fe-Fa)/(ka-ke)
     &            + OaL*(Fe-Fa)*kawz(j)/(ka-ke)
     &            + OaL*ka*(Fe-Fa)*(kewz(j)-kawz(j))/((ka-ke)**2)
     &            + OaL*ka*delta*(kawz(j)*Fa-kewz(j)*Fe)/(ka-ke)
            end do

          end if

          timeL = time(i)
          OaL = Oa
          OL = O
          do j=1,ltau-1
            OawtuL(j)=Oawtau(j)
            OwtauL(j)=Owtau(j)
          end do
          do j=1,M
            OawzL(j)=Oawz(j)
            OwzL(j)=Owz(j)
          end do

          OaVsav(i) = Oa*V
          Osave(i) = O

*         Concentration measurement?
          if (observ(i).eq.0) then

            S=tau(ltau)
            call normalh(conc(i),O,S,2.d0,f,fwO,fwS)
            P = PL*f
            do j=1,ltau-1
              Pwtau(j) = PL*fwO*Owtau(j) + PwtauL(j)*f
            end do
            Pwtau(ltau) = PL*fwS + PwtauL(ltau)*f
            do j=1,M
              Pwz(j) = PL*fwO*Owz(j) + PwzL(j)*f
            end do

            PL = P
            do j=1,ltau
              PwtauL(j) = Pwtau(j)
            end do
            do j=1,M
              PwzL(j) = Pwz(j)
            end do

          end if

        end do

        d = P
        do j=1,ltau
          dwtau(j) = Pwtau(j)
        end do
        do j=1,M
          dwz(j) = Pwz(j)
        end do

        return
      end if

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

      if (isw.eq.3) then
        write(9,9001) id,(z(j),j=1,M)
        return
      end if

*     If isw=4, print plotting information.

      if (isw.eq.4) then

        if (iunit9.ne.1) then
          open (unit=14,file='data.dat',
     &      status='unknown',form='formatted')

          do i=1,mcols
            if (ctsave(i).eq.0) go to 20
            write(14,14002)
     &      i,ctsave(i),idsave(i),
     &      time(i),conc(i),dose(i),
     &      form(i),observ(i),change(i),steady(i),
     &      dinter(i),age(i),height(i),weight(i),
     &      race(i),smoker(i),lush(i),
     &      heart(i),creatn(i),
     &      album(i),glyco(i)
          end do
 20       continue 
          close(14)
        end if
        
        if (iunit9.eq.1) then
      
          fname=outfil//'            '

          j=0
          do i=1,12
            if (fname(i:i).eq.'.') then
              j=i
            end if
            if ((fname(i:i).eq.' ').and.(j.eq.0)) then
              j=i
              fname(j+1:j+1)='f'
            end if
          end do
    
          if (j.eq.0) then
            fname(9:9)='.'
            j=9
          else 
            fname(j:j)='.'
          end if
          
          fname(j+1:j+3)='fit'
          
          if (fname.eq.outfil) then
            fname(j+1:j+3)='ftt'
          end if
          
          open (unit=14,file=fname,status='unknown',form='formatted')
  
          do i=1,mcols
            if (ctsave(i).gt.0.d0) then
              write(14,14003)
     &        i,ctsave(i),idsave(i),
     &        time(i),conc(i),
     &        Osave(i),OaVsav(i),
     &        dose(i),dinter(i)
            end if
          end do

          close(14)
          
          open (unit=14,file='covar.dat',
     &      status='unknown',form='formatted')
     
          do j=1,ncases
            i=point(j)
            write(14,14002)
     &        obs(j),count(j),j,
     &        time(i),conc(i),dose(i),
     &        form(i),observ(i),change(i),steady(i),
     &        dinter(i),age(i),height(i),weight(i),
     &        race(i),smoker(i),lush(i),
     &        heart(i),creatn(i),
     &        album(i),glyco(i)
          end do
          
          close(14)
        end if
        
      return
      end if

      return

3901  format(1x,'error, density, error reading data')
3902  format(1x,'error, density, ',12a,', ',32a)
9001  format(i5,100e30.18)
14001 format(
     &  i4,
     &  f8.2,f4.1,f5.0,
     &  4i2,
     &  f3.0,f3.0,f3.0,f4.0,
     &  3i2,
     &  2i2,
     &  f4.1,f5.0
     &  )
14002 format(
     &  1x,i5,i5,i4,
     &  f8.2,f4.1,f5.0,
     &  4i2,
     &  f4.0,f4.0,f4.0,f5.0,
     &  3i2,
     &  2i2,
     &  f4.1,f5.0
     &  )
14003  format(
     &  1x,i5,i5,i4,
     &  f9.2,f5.1,
     &  f11.5,f13.5,
     &  f6.0,f5.0
     &  )
      end
