C--------------------------------------------------------------------------
C PROGRAM:    a5055.for
C DYNAMIC MODEL: 
C          dT/dt  = lambda-rho*T-[1-gamma(t)]*k*T*V
C          dTp/dt = [1-gamma(t)]*k*T*V - delta*Tp
C          dV/dt  = N*delta*Tp - c*V
C FORM/INPUT(S): hyper-parameters: alpha,beta,eta,big_GAMMA, big_OMAGA,nu
C                C12h, IC50(t), A(t) gamma(t),
C                Data--Viral loads, IC50, PK and adherence
C OUTPUT(S):  estimate parameters: phi,c, delta,lambda,rho,N,k for both po[ulation and
C                                  individuals
C PROTOCOL:  1. THE Runge-Kutta-Verner METHOD OF DIVPRK IS USED
C               TO SOLVE DIFFERENTIAL EQUATIONS
C            2. TO ESTIMATE PARAMETERS USING BAYESIAN METHOD VIA MCMC
C Original program:   07/05/02 Y.Huang
C Modification history: 07/05/05 Y.Huang
C----------------------------------------------------------------------------
C    ****************************************************
C    *(1) MAIN PROGRAM
C    ****************************************************
      program A5055
      implicit double precision (a-h, o-z)
      integer i,j,k,n_sub,n_obs,np,lda,ldainv,nra,nca
      integer nx,ny,ipath,nr,iseed
      parameter(n_sub=42,n_obs=11,n_par=7, nn=328,
     +          n_tot=1+n_par+(n_par*(1+n_par)/2))
      parameter(n_obs10=10,n_obs9=9,n_obs8=8)
      parameter(n_obs7=7,n_obs6=6,n_obs5=5,n_obs4=4,n_obs3=3)
      parameter(nmcmc=100000,nsmcmc=20000,nk=4,
     +          length=(nmcmc-nsmcmc)/nk)
      dimension m(n_sub)
      dimension x(n_sub, n_obs),z(n_sub,n_obs),tp(n_obs+1,n_sub)
      dimension a_mu(n_par),smatinv(n_par,n_par),theta(n_sub,n_par)
      dimension eta(n_par),a_laminv(n_par,n_par),omeginv(n_par,n_par)
      dimension c(n_par),cc(n_par),cct(n_par),tvc(n_par)
      dimension d(n_par,n_par),e(n_par,n_par),estexp(n_par)
      dimension a_mulower(n_par),a_muupper(n_par),a_muexp(n_par)
      dimension a_mumedian(n_par),sd(n_par),sdlog(n_par),vinlog(n_par)
      dimension t_mu(length,n_par),tmu1(length),tmu2(length)
      dimension dinv(n_par,n_par),einv(n_par,n_par), AD(2)
      dimension thetap(n_par),propcov(n_par,n_par),ss(n_sub)
      dimension tmp(n_tot),tmpexp(n_par),exptmp(n_tot),est(n_tot)
      dimension Aidv_Cmin(n_sub),Aidv_IC1(n_sub),Aidv_IC2(n_sub),
     +          Ridv(n_sub)
      dimension Artv_Cmin(n_sub),Artv_IC1(n_sub),Artv_IC2(n_sub),
     +          Rrtv(n_sub)

      dimension theta_sum(n_sub,n_par)
      dimension thetaexp(n_sub,n_par)
      dimension theta_exp(n_sub,n_par)
      dimension est_exp(n_sub,n_par)
      dimension t_theta(n_sub,length,n_par)
      dimension atheta(n_sub,length)
      dimension a_theta(length), b_theta(length)
      dimension sd_s(n_sub,n_par)
      dimension sdlog_s(n_sub,n_par)
      dimension vinlog_s(n_sub,n_par)
      dimension a_me(n_sub,n_par)
      dimension a_l(n_sub,n_par)
      dimension a_u(n_sub,n_par)

      EXTERNAL RNUND,DSCAL,DADD,DRNNOA
      external gamma,wish3,mvnor3,dlinrg,dmurrv,dblinf,drnun,rnset,
     +     dsvrgn,drnunf
C                         The values of hyper-parameters
      parameter (alpha=4.5,beta=9.0,a_nu=8.0)
      data eta/1.5,1.1,-1.0,4.6,-2.5,6.9,-11.0/
      data a_laminv/0.001,  0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
     +              0.0,  400, 0.0, 0.0, 0.0, 0.0, 0.0,
     +              0.0,  0.0, 400, 0.0, 0.0, 0.0, 0.0,
     +              0.0,  0.0, 0.0, 400, 0.0, 0.0, 0.0,
     +              0.0,  0.0, 0.0, 0.0, 400, 0.0, 0.0,
     +              0.0,  0.0, 0.0, 0.0, 0.0, 400, 0.0,
     +              0.0,  0.0, 0.0, 0.0, 0.0, 0.0, 1000/
      data omeginv/0.8, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
     +             0.0, 0.4, 0.0, 0.0, 0.0, 0.0, 0.0,
     +             0.0, 0.0, 0.4, 0.0, 0.0, 0.0, 0.0,
     +             0.0, 0.0, 0.0, 0.5, 0.0, 0.0, 0.0,
     +             0.0, 0.0, 0.0, 0.0, 0.1, 0.0, 0.0, 
     +             0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 0.0,
     +             0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.5/

C                         Data
      data m/n_obs8,n_obs9,n_obs4,n_obs9,n_obs6,n_obs8,n_obs3,
     +       n_obs8,n_obs9,n_obs8,n_obs9,n_obs9,n_obs9,n_obs9,
     +       n_obs9,n_obs9,n_obs6,n_obs8,n_obs, n_obs10,n_obs10,
     +       n_obs5,n_obs8,n_obs9,n_obs7,n_obs9,n_obs6,n_obs9,
     +       n_obs6,n_obs3,n_obs6,n_obs8,n_obs9,n_obs7,n_obs9,
     +       n_obs7,n_obs8,n_obs10,n_obs5,n_obs10,n_obs6,n_obs10/
    
      data propcov/0.01, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
     +             0.0, 0.01, 0.0, 0.0, 0.0, 0.0, 0.0,
     +             0.0, 0.0,  0.01,0.0, 0.0, 0.0, 0.0,
     +             0.0, 0.0,  0.0, 0.01,0.0, 0.0, 0.0,
     +             0.0, 0.0,  0.0, 0.0, 0.01,0.0, 0.0,
     +             0.0, 0.0,  0.0, 0.0, 0.0, 0.01,0.0,
     +             0.0, 0.0,  0.0, 0.0, 0.0, 0.0, 0.01/
C
      PARAMETER (MXPARM=50, N=3)
      PARAMETER(MABSE=1,MNS=3000000) 
      DIMENSION PARAM(MXPARM), Y(N)
      DIMENSION V(n_sub, n_obs+1), A_V0(n_sub),A_PHI(n_sub),A_C0(n_sub)
      DIMENSION A_NUM(n_sub),A_DELTA(n_sub),A_LAMBDA(n_sub),A_AK(n_sub)
      DIMENSION A_RHO(n_sub)
      EXTERNAL DIVPRK,SSET,UMACH,FCN

      COMMON /COM/PHI,C0,DELTA,ALAMBDA,RHO,ANUM,AK
      COMMON /COM1/Cmin_idv, C1_50idv, C2_50idv, AR50idv
      COMMON /COM2/Cmin_rtv, C1_50rtv, C2_50rtv, AR50rtv
      COMMON /COM3/ad_idv, ad_rtv
C                           PK and IC50 data for two drugs
C     --- Ridv is individual failure time point
      data Ridv/  105.0, 0.0, 0.0,   88.0, 0.0, 56.0, 0.0,
     +            0.0,   0.0, 105.0, 0.0,  56.0,0.0,  0.0,
     +            0.0,   0.0, 0.0,   0.0,  0.0, 77.0, 56.0, 0.0,
     +            28.0,  0.0, 78.0,  0.0,  28.0, 0.0, 0.0,
     +            71.0,  0.0, 28.0,  0.0,  0.0,  0.0, 86.0,
     +            0.0,   0.0,  54.0, 0.0,  0.0,  0.0/
     
C     --- Aidv_Cmin is c12h data of IDV drug---
      data Aidv_Cmin/1509.0,1828.4,539.9,529.1,80.4,614.0,770.7,
     +               658.9,633.8,2674.5,422.2,828.9,871.9,105.5,
     +      601.2,1775.2,1037.3,1909.6,744.6,2918.8,648.3,305.4,
     +      655.5,566.8,1451.4,113.5,818.5,343.7,1497.1,548.7,
     +      664.7,164.8, 854.3,1561.0,144.3,793.3,221.6,
     +      739.8,679.1,218.5,620.8,332.9/

C     --- IC50 change at baseline and failure time for IDV drug---
      data Aidv_IC1/2.4,9.6,7.9,4.2,2.5,2.8,6.0,2.8,7.5,2.9,4.3,
     +              2.3,9.9,3.3,7.9,2.5,7.1,11.8,3.1,2.6,2.0,2.8,
     +              2.58,3.31,2.76,1.96,3.44,9.63,3.87,4.42,18.53,5.46,
     +              4.91,3.31,8.71,4.97,5.95,10.12,2.88,10.0,2.45,8.10/
     
      data Aidv_IC2/2.4,9.6,7.9,5.6,2.5,3.3,6.0,2.8,7.5,3.1,4.3,
     +              1.3,9.9,3.3,7.9,2.5,7.1,11.8,3.1,1.7,2.4,2.8,
     +              3.68,3.31,2.58,1.96,1.90,
     +              9.63,3.87,4.05,18.53,3.80,
     +              4.91,3.31,8.71,11.96,5.95,
     +              10.12,1.66,10.0,2.45,8.10/
C     --- Ridv is individual failure time point     
      data Rrtv/105.0, 0.0, 0.0,   88.0, 0.0, 56.0, 0.0, 
     +            0.0,   0.0, 105.0, 0.0,  56.0,0.0,  0.0,
     +            0.0,   0.0, 0.0,   0.0,  0.0, 77.0, 56.0, 0.0,
     +            28.0,  0.0, 78.0,  0.0,  28.0, 0.0, 0.0,
     +            71.0,  0.0, 28.0,  0.0,  0.0,  0.0, 86.0,
     +            0.0,   0.0,  54.0, 0.0,  0.0,  0.0/

C     --- Artv_Cmin is c12h data of RTV drug---
      data Artv_Cmin/1160.3,2964.9,854.3,806.2,666.8,399.4,670.0,
     +               649.2,841.8,2515.9,761.4,2124.7,2714.1,581.3,
     +        448.5,4627.3,1462.7,1067.0,499.5,1727.2,971.2,189.5,
     +               3264.2,2511.1,6107.5,968.6,2539.6,
     +               785.9,6043.6,2322.7,2359.3,655.4,
     +               3681.5,5521.0,258.4,2978.6,556.0,
     +               2616.0,3188.6,2547.0,3017.7,1491.9/

C     --- IC50 change at baseline and failure time for RTV drug---
      data Artv_IC1/9.8,13.9,11.8,17.7,9.6,11.2,28.2,11.8,16.0,11.9,7.5,
     +              8.0,8.7,7.7,37.7,97.1,44.5,222.4,7.7,20.0,3.5,14.6,
     +              8.85,6.69,5.04,4.82,12.09,25.18,13.74,6.91,75.68,
     +   13.02,6.26,7.70,26.62,22.66,20.58,35.54,5.68,49.06,7.70,17.19/
     
      data Artv_IC2/9.8,13.9,11.8,37.5,9.6,15.8,28.2,11.8,16.0,14.7,7.5,
     +              5.8,8.7,7.7,37.7,97.1,44.5,222.4,7.7,10.0,7.1,14.6,
     +         12.81,6.69,9.28,4.82,6.12,25.18,13.74,13.60,75.68,16.26,
     +         6.26,7.70,26.62,55.04,20.58,35.54,5.54,49.06,7.70,17.19/

      open (unit=9,file="a5055.out")
      open (unit=4,file="a5055.dat")
      open (unit=5,file="a5055.par.dat")
      open (unit=6,file="a5055.mu.dat")
      open (unit=7,file="a5055.sample.dat")

      open (unit=2,file="V42.dat")
      open (unit=3,file="day42.dat")

C                                 Read rna and time data
      DO 2 i=1,n_sub
           READ(2,*) (x(i,j),j=1,n_obs)
2     CONTINUE
      DO 3 i=1,n_obs+1
           READ(3,*) (tp(i,j),j=1,n_sub)
3     CONTINUE

      ISEED=12344
      CALL RNSET(ISEED)
      CALL UMACH(2,NOUT)


c        Setup firt value
      do 4 i=1,n_tot
         tmp(i)=0.0
         tmpexp(i)=0.0
4     continue

      do 41 i=1,n_par  
        do 42 l=1,n_sub
           theta_sum(l,i)=0.0
           thetaexp(l,i)=0.0
42      continue
41    continue 

      do 7 i=1,n_sub
         do 6 j=1,m(i)
            z(i,j)=log10(x(i,j))
6       continue
7     continue

C                                 Specifications for optiona in ODE
      TOL=0.01
      CALL SSET(MXPARM,0.0,PARAM,1)
      PARAM(10)=MABSE
      PARAM(4)=MNS  

C                                  initialize
      DO 61 k=1,n_sub
         A_V0(k)=x(k,1)
61    CONTINUE
c      data A_PHI/1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,
c     +           1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1/
c      data A_PHI/12,12,12,12,12,12,12,12,12,12,12,12,12,12,12,
c     +           12,12,12,12,12,12,12,12,12,12,12,12,12,12,12,
c     +           12,12,12,12,12,12,12,12,12,12,12,12/

      data A_PHI/4.5,4.5,4.5,4.5,4.5,4.5,4.5,4.5,4.5,4.5,4.5,
     +           4.5,4.5,4.5,4.5,4.5,4.5,4.5,4.5,4.5,4.5,4.5,
     +           4.5,4.5,4.5,4.5,4.5,4.5,4.5,4.5,4.5,4.5,
     +           4.5,4.5,4.5,4.5,4.5,4.5,4.5,4.5,4.5,4.5/
      data A_C0/3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,
     +          3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3/
      data A_DELTA/0.37,0.37,0.37,0.37,0.37,0.37,0.37,0.37,
     +             0.37,0.37,0.37,0.37,0.37,0.37,0.37,0.37,
     +             0.37,0.37,0.37,0.37,0.37,3.7,
     +             0.37,0.37,0.37,0.37,0.37,0.37,0.37,0.37,
     +             0.37,0.37,0.37,0.37,0.37,0.37,0.37,0.37,
     +             0.37,0.37,0.37,0.37/
      data A_LAMBDA/100,100,100,100,100,100,100,100,
     +              100,100,100,100,100,100,100,100,
     +              100,100,100,100,100,100,
     +              100,100,100,100,100,100,100,100,100,100,
     +              100,100,100,100,100,100,100,100,100,100/
      data A_RHO/0.08,0.08,0.08,0.08,0.08,0.08,0.08,0.08,
     +           0.08,0.08,0.08,0.08,0.08,0.08,0.08,0.08,
     +           0.08,0.08,0.08,0.08,0.08,0.08,0.08,0.08,
     +           0.08,0.08,0.08,0.08,0.08,0.08,0.08,0.08,
     +           0.08,0.08,0.08,0.08,0.08,0.08,0.08,0.08,
     +           0.08,0.08/
      data A_AK/0.00001,0.00001,0.00001,0.00001,0.00001,
     +          0.00001,0.00001,0.00001,0.00001,0.00001,
     +          0.00001,0.00001,0.00001,0.00001,0.00001,
     +          0.00001,0.00001,0.00001,0.00001,0.00001,
     +          0.00001,0.00001,0.00001,0.00001,0.00001,
     +          0.00001,0.00001,0.00001,0.00001,0.00001, 
     +          0.00001,0.00001,0.00001,0.00001,0.00001, 
     +          0.00001,0.00001,0.00001,0.00001,0.00001,  
     +          0.00001,0.00001/

      do 8 i=1,n_sub
         theta(i,1)=log(A_PHI(i))
         theta(i,2)=log(A_C0(i))
         theta(i,3)=log(A_DELTA(i))
         theta(i,4)=log(A_LAMBDA(i))
         theta(i,5)=log(A_RHO(i))
         theta(i,6)=log((A_V0(i)+A_RHO(i)/A_AK(i))*
     +                  (A_C0(i)/A_LAMBDA(i)))
         theta(i,7)=log(A_AK(i))
8     continue
      siginv=1.0
      do 10 i=1,n_par
         do 9 j=1,n_par
            smatinv(i,j)=0.0
9        continue
10    continue
      smatinv(1,1)=5.0
      smatinv(2,2)=5.0
      smatinv(3,3)=5.0
      smatinv(4,4)=0.5
      smatinv(5,5)=1.0
      smatinv(6,6)=1.0
      smatinv(7,7)=0.5
      a_mu(1)=0.0
      a_mu(2)=1.0
      a_mu(3)=-1.0
      a_mu(4)=6.5
      a_mu(5)=-3.5 
      a_mu(6)=4.5
      a_mu(7)=-12
      acc=0.0
      n1=0
      n2=0
      n3=0
      Vlimit=25.0
C                         Start MCMC
      do 999 imcmc=1,nmcmc                 
      a0=alpha+nn/2
99    continue
      acc=acc-n3
      n3=0
      do 13 i=1,n_sub
         ss(i)=0.0
         PHI=exp(theta(i,1))
         C0=exp(theta(i,2))
         DELTA=exp(theta(i,3))
         ALAMBDA=exp(theta(i,4))
         RHO=exp(theta(i,5))
         AK=exp(theta(i,7))
         ANUM=exp(theta(i,6))
c         print *, 'i and First est par is:', i
c         write(*,'(7E12.5)') PHI,C0,DELTA,ALAMBDA, RHO,ANUM,AK

         Cmin_idv=Aidv_Cmin(i)
         C1_50idv=Aidv_IC1(i)
         C2_50idv=Aidv_IC2(i)
         AR50idv=Ridv(i)

         Cmin_rtv=Artv_Cmin(i)
         C1_50rtv=Artv_IC1(i)
         C2_50rtv=Artv_IC2(i)
         AR50rtv=Rrtv(i)
 
         V0=ALAMBDA*ANUM/C0-RHO/AK
         IF (V0 .LE. 0.0) THEN
            n1=n1+1
         END IF

         T0=C0/(AK*ANUM)
         Tp0=C0*V0/(DELTA*ANUM)

         T=tp(1,i)
         Y(1)=T0
         Y(2)=Tp0
         Y(3)=V0        
         V(i,1)=log10(Y(3))
         istep=1
         ss(i)=ss(i)+(z(i,1)-V(i,1))**2
         IDO=1
12       CONTINUE
         istep=istep+1
         j=istep
         TEND=tp(j,i)

C        ****** CALL SUBROUTINE FOR ADHERENCE*****
         CALL ADHEREN(i, TEND, AD)
         ad_idv=AD(1)
         ad_rtv=AD(2)

C        ****** CALL IMSL SUBROUTINE to solve ODE    
         CALL DIVPRK(IDO,N,FCN,T,TEND,TOL,PARAM,Y)
         IF (Y(3) .LT. Vlimit) THEN
               Y(3)=Vlimit
         END IF
         V(i,j)=log10(Y(3))
         IF (tp(j,i) .LE. tp(m(i),i)) THEN
              IF (tp(j,i) .EQ. tp(m(i),i)) IDO=3 
              ss(i)=ss(i)+(z(i,j)-V(i,j))**2
              GO TO 12
         END IF
13    continue
      sum_ssi=0
      do 133 i=1,n_sub
         sum_ssi=sum_ssi+ss(i)
133   continue
      b=1.0/(1.0/beta+0.5*sum_ssi)
      call gamma(siginv,a0,b)

      do 15 i=1,n_par
        do 14 j=1,n_par
            dinv(i,j)=dble(n_sub)*smatinv(i,j)+a_laminv(i,j)
14      continue
15    continue
      np=n_par
      lda=n_par
      ldainv=n_par
      call dlinrg(np,dinv,lda,d,ldainv)
      nra=n_par
      nca=n_par
      lda=n_par
      nx=n_par
      ny=n_par
      ipath=2
      call dmurrv(nra,nca,a_laminv,lda,nx,eta,ipath,ny,cc)
      do 18 i=1,n_sub
         do 16 j=1,n_par
            tvc(j)=theta(i,j)
16       continue
         call dmurrv(nra,nca,smatinv,lda,nx,tvc,ipath,ny,cct)
         do 17 j=1,n_par
            cc(j)=cc(j)+cct(j)
17       continue
18    continue
      call dmurrv(nra,nca,d,lda,nx,cc,ipath,ny,c)
      call mvnor3(a_mu,c,d)
      f=dble(n_sub)+a_nu
      do 21 i=1,n_par
         do 20 j=1,n_par
            einv(i,j)=omeginv(i,j)
            do 19 k=1,n_sub
             einv(i,j)=einv(i,j)+(theta(k,i)-a_mu(i))*
     +                 (theta(k,j)-a_mu(j))
19          continue
20       continue
21    continue
      call dlinrg(np,einv,lda,e,ldainv)
      call wish3(smatinv,e,f)

      do 510 i=1,n_sub
         do 520 j=1,n_par
            tvc(j)=theta(i,j)
 520     continue
         call mvnor3(thetap,tvc,propcov)
         prob=0.5*siginv*ss(i)
         do 530 j=1,n_par
            tvc(j)=theta(i,j)-a_mu(j)
 530     continue
         prob=prob+0.5*dblinf(nra,nca,smatinv,lda,tvc,tvc)
         do 540 j=1,n_par
            tvc(j)=thetap(j)-a_mu(j)
 540     continue
         prob=prob-0.5*dblinf(nra,nca,smatinv,lda,tvc,tvc)
C
         A_PHI(i)=exp(thetap(1))
         A_C0(i)=exp(thetap(2))
         A_DELTA(i)=exp(thetap(3))
         A_LAMBDA(i)=exp(thetap(4))
         A_RHO(i)=exp(thetap(5))
         A_NUM(i)=exp(thetap(6))
         A_AK(i)=exp(thetap(7))

         PHI=A_PHI(i)
         C0=A_C0(i)
         DELTA=A_DELTA(i)
         ALAMBDA=A_LAMBDA(i)
         RHO=A_RHO(i)
         ANUM=A_NUM(i)
         AK=A_AK(i)
c         print *, 'i and Second est par is', i
c         write(*,'(7E12.5)') PHI,C0,DELTA,ALAMBDA, RHO,ANUM,AK


         Cmin_idv=Aidv_Cmin(i)
         C1_50idv=Aidv_IC1(i)
         C2_50idv=Aidv_IC2(i)
         AR50idv=Ridv(i)
     
         Cmin_rtv=Artv_Cmin(i)
         C1_50rtv=Artv_IC1(i)
         C2_50rtv=Artv_IC2(i)
         AR50rtv=Rrtv(i)

         V0=ALAMBDA*ANUM/C0-RHO/AK
         IF (V0 .LE. 0.0) THEN
C            V0=A_V0(i)
            n2=n2+1
         END IF
         IF (V0 .LE. 1.0) GO TO 99
         T0=C0/(AK*ANUM)
         Tp0=C0*V0/(DELTA*ANUM)

         T=tp(1,i)
         Y(1)=T0
         Y(2)=Tp0
         Y(3)=V0
         V(i,1)=log10(Y(3))
         istep=1
         prob=prob-0.5*siginv*(z(i,1)-V(i,1))**2
         IDO=1        
550      CONTINUE
         istep=istep+1
         j=istep
         TEND=tp(j,i)

C        ****** CALL SUBROUTINE FOR ADHERENCE*****
         CALL ADHEREN(i, TEND, AD)
         ad_idv=AD(1)
         ad_rtv=AD(2)
C        ****** CALL IMSL SUBROUTINE to solve ODE           
         CALL DIVPRK(IDO,N,FCN,T,TEND,TOL,PARAM,Y)
         IF (Y(3) .LT. Vlimit) THEN
               Y(3)=Vlimit
         END IF
         V(i,j)=log10(Y(3))
         IF (tp(j,i) .LE. tp(m(i),i)) THEN
              IF (tp(j,i) .EQ. tp(m(i),i)) IDO=3
              prob=prob-0.5*siginv*(z(i,j)-V(i,j))**2
              GO TO 550
         END IF
         prob=exp(prob)
         nr=1
         call drnun(nr,uni)
         if (uni .le. prob) then
            do 560 j=1,n_par
               theta(i,j)=thetap(j)
 560        continue
            acc=acc+1.0
            n3=n3+1
         endif
 510  continue
           n3=0
      if (imcmc .gt. nsmcmc .and. mod(imcmc,nk) .eq. 0) then
         do 721 i=1,n_par
                tmp(i)=tmp(i)+a_mu(i)
                tmpexp(i)=tmpexp(i)+exp(a_mu(i))
                a_muexp(i)=exp(a_mu(i))
                t_mu((imcmc-nsmcmc)/nk,i)=exp(a_mu(i))
721      continue
         tmp(8)=tmp(8)+smatinv(1,1)
         tmp(9)=tmp(9)+smatinv(1,2)
         tmp(10)=tmp(10)+smatinv(1,3)
         tmp(11)=tmp(11)+smatinv(1,4)
         tmp(12)=tmp(12)+smatinv(1,5)
         tmp(13)=tmp(13)+smatinv(1,6)
         tmp(14)=tmp(14)+smatinv(1,7)
         tmp(15)=tmp(15)+smatinv(2,2)
         tmp(16)=tmp(16)+smatinv(2,3)
         tmp(17)=tmp(17)+smatinv(2,4)
         tmp(18)=tmp(18)+smatinv(2,5)
         tmp(19)=tmp(19)+smatinv(2,6)
         tmp(20)=tmp(20)+smatinv(2,7)
         tmp(21)=tmp(21)+smatinv(3,3)
         tmp(22)=tmp(22)+smatinv(3,4)
         tmp(23)=tmp(23)+smatinv(3,5)
         tmp(24)=tmp(24)+smatinv(3,6)
         tmp(25)=tmp(25)+smatinv(3,7)
         tmp(26)=tmp(26)+smatinv(4,4)
         tmp(27)=tmp(27)+smatinv(4,5)
         tmp(28)=tmp(28)+smatinv(4,6)
         tmp(29)=tmp(29)+smatinv(4,7)
         tmp(30)=tmp(30)+smatinv(5,5)   
         tmp(31)=tmp(31)+smatinv(5,6)
         tmp(32)=tmp(32)+smatinv(5,7)
         tmp(33)=tmp(33)+smatinv(6,6)
         tmp(34)=tmp(34)+smatinv(6,7)
         tmp(35)=tmp(35)+smatinv(7,7)
         tmp(36)=tmp(36)+siginv
        
         write(6, 881) (a_mu(k),k=1,n_par)
         do 741 i=1,n_par
             ndiff=(imcmc-nsmcmc)/nk
             est(i)=tmp(i)/ndiff
741      continue
         write(7, 882) imcmc-nsmcmc, (est(k), k=1,n_par)

C        *******ADD3*******
         do 7211 j=1,n_par
            do 7212 l=1,n_sub
                theta_sum(l,j)=theta_sum(l,j)+theta(l,j)
                thetaexp(l,j)=thetaexp(l,j)+exp(theta(l,j))
                t_theta(l,(imcmc-nsmcmc)/nk,j)=exp(theta(l,j))
7212        continue
7211     continue
      endif
c       if (mod(imcmc,1000) .eq. 0) print *,imcmc
999   continue
      n_lower1=length*25/1000  
      n_lower2=n_lower1+1  
      n_upper1=length*975/1000  
      n_upper2=n_upper1+1  
      do 805 j=1,n_tot
         tmp(j)=tmp(j)/length
805   continue
      do 8050 j=1,n_par
         tmpexp(j)=tmpexp(j)/length
         exptmp(j)=exp(tmp(j))
8050  continue

      do 806 i=1,n_par
         sd(i)=0.0
         sdlog(i)=0.0
         do 807 k=1,length
               tmu1(k)=t_mu(k,i)
               sd(i)=sd(i)+(t_mu(k,i)-tmpexp(i))**2
               sdlog(i)=sdlog(i)+(log(t_mu(k,i))-tmp(i))**2
807      continue

         a_mumedian(i)=t_mu(int(length/2),i)
         sd(i)=sqrt(sd(i)/(length-1))
         vinlog(i)=1/(sdlog(i)/(length-1))
         sdlog(i)=sqrt(sdlog(i)/(length-1))
         call dsvrgn(length,tmu1,tmu2)
         a_mulower(i)=0.5*(tmu2(n_lower1)+tmu2(n_lower2))
         a_muupper(i)=0.5*(tmu2(n_upper1)+tmu2(n_upper2))
         estexp(i)=exp(tmp(i)+0.5*sdlog(i)**2)
806   continue

C     ******ADD4*******
      do 8051 j=1,n_par
         do 8052 l=1,n_sub  
           theta_sum(l,j)=theta_sum(l,j)/length
           theta_exp(l,j)=exp(theta_sum(l,j))
           thetaexp(l,j)=thetaexp(l,j)/length
8052     continue
8051   continue 

      do 8061 j=1,n_par
        do 8062 l=1,n_sub
         sd_s(l,j)=0.0
         sdlog_s(l,j)=0.0

         do 8071 k=1,length
          atheta(l,k)=t_theta(l,k,j)
          sd_s(l,j)=sd_s(l,j)+(atheta(l,k)-thetaexp(l,j))**2
          sdlog_s(l,j)=sdlog_s(l,j)+(log(atheta(l,k))-
     +                 theta_sum(l,j))**2
8071     continue

         vinlog_s(l,j)=1/(sdlog_s(l,j)/(length-1))
         a_me(l,j)=atheta(l,int(length/2))
         sd_s(l,j)=sqrt(sd_s(l,j)/(length-1))
         sdlog_s(l,j)=sqrt(sdlog_s(l,j)/(length-1))
        
      do 8063 k=1,length
         a_theta(k)=atheta(l,k)
8063  continue
         call dsvrgn(length,a_theta,b_theta)
         a_l(l,j)=0.5*(b_theta(n_lower1)+b_theta(n_lower2))
         a_u(l,j)=0.5*(b_theta(n_upper1)+b_theta(n_upper2))

         est_exp(l,j)=exp(theta_sum(l,j)+0.5*sdlog_s(l,j)**2)
8062     continue
8061    continue

      acc=acc/(nmcmc*n_sub)
      write (9,870) nmcmc,nsmcmc,length,nk,ISEED
      write (9,883) acc, n1,n2,ISEED, TOL
      write (9,871) (exptmp(k),k=1,7)
      write (9,8710) (tmpexp(k),k=1,7)
      write (9,8711) estexp
      write (9,872) a_mulower  
      write (9,873) a_muupper
      write (9,874) a_mumedian
      write (9,875) sd
      write (9,8751) sdlog
      write (9,8752) vinlog
      write (9,876) tmp

      do 100 l=1,n_sub
      write (9,884) l
      write (9,871)  (theta_exp(l,j),j=1,n_par)
      write (9,8710) (thetaexp(l,j),j=1,n_par)
      write (9,8711) (est_exp(l,j),j=1,n_par)
      write (9,872)  (a_l(l,j),j=1,n_par)
      write (9,873)  (a_u(l,j),j=1,n_par)
      write (9,874)  (a_me(l,j),j=1,n_par)
      write (9,875)  (sd_s(l,j),j=1,n_par)
      write (9,8751) (sdlog_s(l,j),j=1,n_par)
      write (9,8752) (vinlog_s(l,j),j=1,n_par)
      write (9,876)  (theta_sum(l,j),j=1,n_par)

      write (4,770) (est_exp(l,j),j=1,n_par)
      write (5,771) (theta_exp(l,j),j=1,n_par)
100   continue
      close (unit=9)
      close (unit=2)
      close (unit=3)
      close (unit=4)      
      close (unit=5)
      close (unit=6)
      close (unit=7)
870   format (1X, i8, 1X,i6,1X,i8,1X,i3,1X,i7)
871   format('exptmp=', 5f10.3,2x,f10.3,1x,f11.7)
8710  format('tmpexp=', 5f10.3,2x,f10.3,1x,f11.7)
8711  format('estexp=', 5f10.3,2x,f10.3,1x,f11.7)
872   format('mulow=',  1x,5f10.3,2x,f10.3,1x,f11.7)
873   format('muupp=',  1x,5f10.3,2x,f10.3,1x,f11.7)
874   format('median=', 5f10.3,2x,f10.3,1x,f11.7)
875   format('sd=',     4x,5f10.3,2x,f10.3,1x,f11.7)
8751  format('sdlog=',  1x,5f10.3,2x,f10.3,1x,f11.7)
8752  format('vinlog=', 5f10.3,2x,f10.3,1x,f11.3)
876   format('tmp=',    3x,5f10.3,2x,f10.3,1x,f11.4)

770   format(5f10.3,2x,f10.3,1x,f11.7)
771   format(5f10.3,2x,f10.3,1x,f11.7)
772   format(5f10.3,2x,f10.3,1x,f11.7)
773   format(5f10.3,2x,f10.3,1x,f11.3)
774   format(5f10.3,2x,f10.3,1x,f11.6)
881   format (7(1X,F10.4))
882   format(i6,1x,7(1X,F10.4)) 
883   format(f7.4,1X,i6,1X,i6,1X,i7,1x, F8.5)
884   format('nosub=', i3)
990   format('1st=', f10.2,1x,f10.4)
991   format('2nd=', f10.2,1x,f10.4,f10.5,1x,f14.9)

      stop
      end
C   ****************************************************
C   * SUBROUTINES AND FUNCTIONS
C   ****************************************************
C (1)                                To generate gamma samples
      subroutine gamma (x,alpha,beta)
      implicit double precision (a-h, o-z)
      external drngam
      nr=1
      call drngam(nr,alpha,x)
      x=x*beta
      return
      end
C (2)                                To generate Wishart samples
      subroutine wish3(v,omega,a_nu)
      implicit double precision (a-h, o-z)
      integer n,ldr,irank,nr,nra,nca,lda,nrb,ncb,ldb,nrc,ncc,ldc,nb
      parameter(n_par=7)
      dimension v(n_par,n_par),omega(n_par,n_par)
      dimension chol(n_par,n_par),awish(n_par,n_par),hwish(n_par,n_par)
      intrinsic sqrt
      external drnnor,drngam,dchfac,dmrrrr,dmxtxf,dmach
      nr=1
      n=n_par
      lda=n_par
      ldr=n_par
      tol=100.0*dmach(4)
      call dchfac(n,omega,lda,tol,irank,chol,ldr)
      call drnnor(nr,awish(1,2))
      call drnnor(nr,awish(1,3))
      call drnnor(nr,awish(1,4))
      call drnnor(nr,awish(1,5))
      call drnnor(nr,awish(1,6))
      call drnnor(nr,awish(1,7))
      call drnnor(nr,awish(2,3))
      call drnnor(nr,awish(2,4))
      call drnnor(nr,awish(2,5))
      call drnnor(nr,awish(2,6))
      call drnnor(nr,awish(2,7))
      call drnnor(nr,awish(3,4))
      call drnnor(nr,awish(3,5))
      call drnnor(nr,awish(3,6))
      call drnnor(nr,awish(3,7))
      call drnnor(nr,awish(4,5))
      call drnnor(nr,awish(4,6))
      call drnnor(nr,awish(4,7))
      call drnnor(nr,awish(5,6))
      call drnnor(nr,awish(5,7))
      call drnnor(nr,awish(6,7))  
      shape=a_nu/2.0
      call drngam(nr,shape,awish(1,1))
      awish(1,1)=sqrt(2.0*awish(1,1))
      shape=(a_nu-1.0)/2.0
      call drngam(nr,shape,awish(2,2))
      awish(2,2)=sqrt(2.0*awish(2,2))
      shape=(a_nu-2.0)/2.0
      call drngam(nr,shape,awish(3,3))
      awish(3,3)=sqrt(2.0*awish(3,3))
      shape=(a_nu-3.0)/2.0
      call drngam(nr,shape,awish(4,4))
      awish(4,4)=sqrt(2.0*awish(4,4))
      shape=(a_nu-4.0)/2.0
      call drngam(nr,shape,awish(5,5))
      awish(5,5)=sqrt(2.0*awish(5,5))
      shape=(a_nu-5.0)/2.0
      call drngam(nr,shape,awish(6,6))
      awish(6,6)=sqrt(2.0*awish(6,6))
      shape=(a_nu-6.0)/2.0
      call drngam(nr,shape,awish(7,7))
      awish(7,7)=sqrt(2.0*awish(7,7))
      awish(2,1)=0.0
      awish(3,1)=0.0
      awish(3,2)=0.0
      awish(4,1)=0.0
      awish(4,2)=0.0
      awish(4,3)=0.0
      awish(5,1)=0.0
      awish(5,2)=0.0
      awish(5,3)=0.0
      awish(5,4)=0.0
      awish(6,1)=0.0
      awish(6,2)=0.0
      awish(6,3)=0.0
      awish(6,4)=0.0
      awish(6,5)=0.0
      awish(7,1)=0.0
      awish(7,2)=0.0
      awish(7,3)=0.0
      awish(7,4)=0.0
      awish(7,5)=0.0
      awish(7,6)=0.0
      nra=n_par
      nca=n_par
      lda=n_par
      nrb=n_par
      ncb=n_par
      ldb=n_par
      nrc=n_par
      ncc=n_par
      ldc=n_par
      call dmrrrr(nra,nca,awish,lda,nrb,ncb,chol,ldb,nrc,ncc,hwish,ldc)
      nb=n_par
      call dmxtxf(nra,nca,hwish,lda,nb,v,ldb)
      return
      end
C (3)                              To generate multivariate normal samples
      subroutine mvnor3(x,a_mu,sigma)
      implicit double precision (a-h, o-z)
      integer i,irank,k,ldr,ldrsig,nr,lda
      parameter(n_par=7)
      dimension  x(n_par),a_mu(n_par),sigma(n_par,n_par),
     +           rsig(n_par,n_par)
      external dchfac,drnmvn,dmach
      k=n_par
      lda=n_par
      ldrsig=n_par
      tol=100.0*dmach(4)
      call dchfac(k,sigma,lda,tol,irank,rsig,ldrsig)
      nr=1
      ldr=1
      call drnmvn(nr,k,rsig,ldrsig,x,ldr)
      do 10 i=1,n_par
         x(i)=x(i)+a_mu(i)
 10   continue
      return
      end
C                                 To solve differential equations
C (4)
      SUBROUTINE FCN(N, T, Y, YPRIME)
      IMPLICIT DOUBLE PRECISION (A-H, O-Z)
      INTEGER N
      DIMENSION Y(N), YPRIME(N)
      COMMON /COM/PHI,C0,DELTA,ALAMBDA,RHO,ANUM,AK
      COMMON /COM1/Cmin_idv, C1_50idv, C2_50idv, AR50idv
      COMMON /COM2/Cmin_rtv, C1_50rtv, C2_50rtv, AR50rtv
      COMMON /COM3/ad_idv, ad_rtv
      C50idv=C50FUN(AR50idv,C1_50idv,C2_50idv,T)
      C50rtv=C50FUN(AR50rtv,C1_50rtv,C2_50rtv,T)
      adh1=ad_idv
      adh2=ad_rtv
      DNOM=Cmin_idv*adh1/C50idv+Cmin_rtv*adh2/C50rtv
      DEN=PHI+DNOM
      IF (T .EQ. 0.0) THEN
          GAMMA=0.0
      ELSE 
          GAMMA=DNOM/DEN
      END IF
      YPRIME(1)=ALAMBDA-RHO*Y(1)-(1-GAMMA)*AK*Y(1)*Y(3)
      YPRIME(2)=(1-GAMMA)*AK*Y(1)*Y(3)-DELTA*Y(2)
      YPRIME(3)=ANUM*DELTA*Y(2)-C0*Y(3)
      RETURN
      END

C (5) Adherence
      SUBROUTINE ADHEREN(i, TEND, AD)
      IMPLICIT DOUBLE PRECISION (A-H, O-Z)
      DIMENSION AD(2)
C         **(1)**
         IF (i .EQ. 1) THEN 
            IF (TEND .LE. 56.0) THEN
                ad_idv=0.79
                ad_rtv=0.77  
            ELSE IF (TEND .LE. 149.0) THEN
                ad_idv=0.65 
                ad_rtv=0.54
            ELSE IF (TEND .LE. 169.0) THEN
                ad_idv=0.83
                ad_rtv=1.00
            ELSE
                ad_idv=1.00
                ad_rtv=1.00
            END IF
C         **(2)** 
         ELSE IF (i .EQ. 2) THEN
            IF (TEND .LE. 13.0) THEN
                ad_idv=1.00
                ad_rtv=1.00
            ELSE IF (TEND .LE. 26.0) THEN
                ad_idv=0.88
                ad_rtv=0.77
            ELSE IF (TEND .LE. 111.0) THEN
                ad_idv=1.00
                ad_rtv=1.00
            ELSE IF (TEND .LE. 167.0) THEN
                ad_idv=0.98
                ad_rtv=1.00
            ELSE
                ad_idv=1.00
                ad_rtv=1.00
            END IF
C         **(3)**
         ELSE  IF (i .EQ. 3) THEN
            ad_idv=1.00
            ad_rtv=1.00
C         **(4)**
         ELSE IF (i .EQ. 4) THEN
            ad_idv=1.00
            ad_rtv=1.00
C         **(5)**
         ELSE IF (i .EQ. 5) THEN
            IF (TEND .LE. 34.0) THEN
                ad_idv=0.91
                ad_rtv=0.93
            ELSE
                ad_idv=1.00
                ad_rtv=1.00
            END IF
C         **(6)**
         ELSE IF (i .EQ. 6) THEN
            IF (TEND .LE. 140.0) THEN
                ad_idv=0.96
                ad_rtv=0.85
            ELSE
                ad_idv=1.00
                ad_rtv=1.00
            END IF
C         **(7)**
         ELSE IF (i .EQ. 7) THEN
             ad_idv=1.00
             ad_rtv=1.00
C         **(8)**
         ELSE IF (i .EQ. 8) THEN
            IF (TEND .LE. 64.0) THEN
                ad_idv=0.95
                ad_rtv=0.95
            ELSE
                ad_idv=1.00
                ad_rtv=1.00
            END IF
C         **(9)**
         ELSE IF (i .EQ. 9) THEN
            IF (TEND .LE. 29.0) THEN
                ad_idv=1.00
                ad_rtv=1.00
            ELSE IF (TEND .LE. 56.0) THEN
                ad_idv=1.00
                ad_rtv=0.89
            ELSE IF (TEND .LE. 86.0) THEN
                ad_idv=1.00
                ad_rtv=1.00
            ELSE IF (TEND .LE. 114.0) THEN
                ad_idv=1.00
                ad_rtv=0.89
            ELSE 
                ad_idv=1.00
                ad_rtv=1.00
            END IF
C         **(10)**
         ELSE IF (i .EQ. 10) THEN
             ad_idv=1.00
             ad_rtv=1.00
C         **(11)**
         ELSE IF (i .EQ. 11) THEN
            IF (TEND .LE. 14.0) THEN
                ad_idv=0.97
                ad_rtv=1.00
            ELSE 
                ad_idv=1.00
                ad_rtv=1.00
            END IF
C         **(12)**
         ELSE IF (i .EQ. 12) THEN
             ad_idv=1.00
             ad_rtv=1.00  
C         **(13)**
         ELSE IF (i .EQ. 13) THEN
            IF (TEND .LE. 112.0) THEN
                ad_idv=1.00
                ad_rtv=1.00
            ELSE IF (TEND .LE. 141.0) THEN
                ad_idv=0.98
                ad_rtv=0.98
            ELSE
                ad_idv=1.00
                ad_rtv=1.00
            END IF
C         **(14)**
         ELSE IF (i .EQ. 14) THEN
            IF (TEND .LE. 57.0) THEN
                ad_idv=1.00
                ad_rtv=1.00
            ELSE IF (TEND .LE. 108.0) THEN
                ad_idv=1.00
                ad_rtv=0.94
            ELSE IF (TEND .LE. 128.0) THEN
                ad_idv=1.00
                ad_rtv=1.00
            ELSE IF (TEND .LE. 163.0) THEN
                ad_idv=0.99
                ad_rtv=0.99
            ELSE
                ad_idv=1.00
                ad_rtv=1.00
            END IF
C         **(15)**
         ELSE IF (i .EQ. 15) THEN
            IF (TEND .LE. 112.0) THEN
                ad_idv=1.00
                ad_rtv=1.00
            ELSE IF (TEND .LE. 144.0) THEN
                ad_idv=0.98
                ad_rtv=0.98
            ELSE
                ad_idv=1.00
                ad_rtv=1.00
            END IF
C         **(16)**
         ELSE IF (i .EQ. 16) THEN
            IF (TEND .LE. 28.0) THEN
                ad_idv=0.82
                ad_rtv=0.84
            ELSE
                ad_idv=1.00
                ad_rtv=1.00
            END IF
C         **(17)**
         ELSE IF (i .EQ. 17) THEN
            IF (TEND .LE. 34.0) THEN
                ad_idv=1.00
                ad_rtv=1.00
            ELSE IF (TEND .LE. 63.0) THEN
                ad_idv=0.98
                ad_rtv=0.98
            ELSE IF (TEND .LE. 91.0) THEN
                ad_idv=0.96
                ad_rtv=0.91
            ELSE
                ad_idv=1.00
                ad_rtv=1.00
            END IF
C         **(18)**
         ELSE IF (i .EQ. 18) THEN
             ad_idv=1.00  
             ad_rtv=1.00
C         **(19)**
         ELSE IF (i .EQ. 19) THEN
            IF (TEND .LE. 13.0) THEN
                ad_idv=1.00
                ad_rtv=1.00
            ELSE IF (TEND .LE. 28.0) THEN
                ad_idv=0.33
                ad_rtv=0.33
            ELSE 
                ad_idv=1.00
                ad_rtv=1.00
            END IF
C         **(20)**
         ELSE IF (i .EQ. 20) THEN
             ad_idv=1.00
             ad_rtv=1.00  
C         **(21)**
         ELSE IF (i .EQ. 21) THEN
             ad_idv=1.00  
             ad_rtv=1.00  
c         **(22)**
         ELSE IF (i .EQ. 22) THEN
             ad_idv=1.00  
             ad_rtv=1.00 
C         **(23)**
         ELSE IF (i .EQ. 23) THEN
             ad_idv=1.00  
             ad_rtv=1.00
C         **(24)**
         ELSE IF (i .EQ. 24) THEN
             ad_idv=1.00
             ad_rtv=1.00  
C       **(25)**
         ELSE IF (i .EQ. 25) THEN
            IF (TEND .LE. 13.0) THEN
                ad_idv=1.00
                ad_rtv=1.00
            ELSE IF (TEND .LE. 56.0) THEN
                ad_idv=0.99
                ad_rtv=0.99
            ELSE 
                ad_idv=1.00
                ad_rtv=1.00
           END IF
C         **(26)**
         ELSE IF (i .EQ. 26) THEN
            IF (TEND .LE. 14.0) THEN
                ad_idv=0.93
                ad_rtv=1.00
            ELSE IF (TEND .LE. 28.0) THEN
                ad_idv=0.96
                ad_rtv=0.96
            ELSE IF (TEND .LE. 54.0) THEN
                ad_idv=1.00
                ad_rtv=0.98
            ELSE IF (TEND .LE. 82.0) THEN
                ad_idv=0.98
                ad_rtv=0.99
            ELSE IF (TEND .LE. 112.0) THEN
                ad_idv=0.83
                ad_rtv=0.83
            ELSE IF (TEND .LE. 146.0) THEN
                ad_idv=0.72
                ad_rtv=0.68
            ELSE
                ad_idv=1.00
                ad_rtv=1.00
            END IF
C         **(27)**
         ELSE IF (i .EQ. 27) THEN
            IF (TEND .LE. 14.0) THEN
                ad_idv=0.34
                ad_rtv=0.66
            ELSE IF (TEND .LE. 28.0) THEN
                ad_idv=0.75
                ad_rtv=0.77
            ELSE
                ad_idv=1.00
                ad_rtv=1.00
            END IF
C         **(28)**
         ELSE IF (i .EQ. 28) THEN
            IF (TEND .LE. 82.0) THEN
                ad_idv=0.87
                ad_rtv=0.96
            ELSE IF (TEND .LE. 112.0) THEN
                ad_idv=0.93
                ad_rtv=0.95
            ELSE
                ad_idv=1.00
                ad_rtv=1.00
            END IF
C         **(29)**
         ELSE IF (i .EQ. 29) THEN
            IF (TEND .LE. 20.0) THEN
                ad_idv=1.00
                ad_rtv=1.00
            ELSE IF (TEND .LE. 99.0) THEN
                ad_idv=0.97
                ad_rtv=0.99
            ELSE
                ad_idv=1.00
                ad_rtv=1.00
            END IF
C         **(30)**
         ELSE IF (i .EQ. 30) THEN
             ad_idv=1.00  
             ad_rtv=1.00
C         **(31)**
         ELSE IF (i .EQ. 31) THEN
             ad_idv=1.00
             ad_rtv=1.00  
C         **(32)**
         ELSE IF (i .EQ. 32) THEN
             ad_idv=1.00  
             ad_rtv=1.00  
C         **(33)**
         ELSE IF (i .EQ. 33) THEN
            IF (TEND .LE. 112.0) THEN
                ad_idv=1.00
                ad_rtv=1.00
            ELSE IF (TEND .LE. 140.0) THEN
                ad_idv=0.98
                ad_rtv=0.96
            ELSE
                ad_idv=1.00
                ad_rtv=1.00
            END IF
C         **(34)**
         ELSE IF (i .EQ. 34) THEN
            IF (TEND .LE. 20.0) THEN
                ad_idv=1.00
                ad_rtv=1.00
            ELSE IF (TEND .LE. 61.0) THEN
                ad_idv=0.96
                ad_rtv=0.96
            ELSE IF (TEND .LE. 90.0) THEN
                ad_idv=1.00
                ad_rtv=1.00
            ELSE IF (TEND .LE. 119.0) THEN
                ad_idv=0.98
                ad_rtv=0.98
            ELSE IF (TEND .LE. 166.0) THEN
                ad_idv=0.98
                ad_rtv=0.96
            ELSE
                ad_idv=1.00
                ad_rtv=1.00
            END IF
C         **(35)**
         ELSE IF (i .EQ. 35) THEN
            IF (TEND .LE. 26.0) THEN
                ad_idv=1.00
                ad_rtv=1.00
            ELSE IF (TEND .LE. 53.0) THEN
                ad_idv=0.99
                ad_rtv=1.00
            ELSE
                ad_idv=1.00
                ad_rtv=1.00
            END IF
C         **(36)**
         ELSE IF (i .EQ. 36) THEN
             ad_idv=1.00  
             ad_rtv=1.00
C         **(37)**
         ELSE IF (i .EQ. 37) THEN
             ad_idv=1.00
             ad_rtv=1.00  
C         **(38)**
         ELSE IF (i .EQ. 38) THEN
            IF (TEND .LE. 56.0) THEN
                ad_idv=1.00
                ad_rtv=1.00
            ELSE IF (TEND .LE. 84.0) THEN
                ad_idv=0.93
                ad_rtv=0.94
            ELSE IF (TEND .LE. 111.0) THEN
                ad_idv=1.00
                ad_rtv=1.00
            ELSE IF (TEND .LE. 140.0) THEN
                ad_idv=0.98
                ad_rtv=0.98
            ELSE
                ad_idv=1.00
                ad_rtv=1.00
            END IF
C         **(39)**
         ELSE IF (i .EQ. 39) THEN
             ad_idv=1.00
             ad_rtv=1.00
C         **(40)**
         ELSE IF (i .EQ. 40) THEN
             ad_idv=1.00
             ad_rtv=1.00
C         **(41)**
         ELSE IF (i .EQ. 41) THEN
             ad_idv=1.00
             ad_rtv=1.00
C         **(42)**
         ELSE IF (i .EQ. 42) THEN
            IF (TEND .LE. 15.0) THEN
                ad_idv=0.97
                ad_rtv=0.90
            ELSE IF (TEND .LE. 25.0) THEN
                ad_idv=1.00
                ad_rtv=1.00
            ELSE IF (TEND .LE. 56.0) THEN
                ad_idv=0.97
                ad_rtv=0.97
            ELSE IF (TEND .LE. 84.0) THEN
                ad_idv=1.00
                ad_rtv=1.00
            ELSE IF (TEND .LE. 140.0) THEN
                ad_idv=1.00
                ad_rtv=0.99
            ELSE
                ad_idv=1.00
                ad_rtv=1.00
            END IF
         END IF
         AD(1)=ad_idv
         AD(2)=ad_rtv  
      RETURN
      END
C (6)     Function IC50(t)
      FUNCTION C50FUN(D,C1,C2,T)
      IMPLICIT DOUBLE PRECISION (A-H, O-Z)
      IF (D .EQ. 0.0) THEN
         C50=C1
      ELSE IF (T .LE. D) THEN
         C50=C1+T*(C2-C1)/D
      ELSE
         C50=C2
      END IF
      C50FUN=C50
      RETURN  
      END          