A5055 Long-Term Viral Dynamic Data (including PK, Adherence, and Drug Resistance): Hierarchical Bayesian (NLME) Models

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