1! disturbance smoothing algorithm for simulation
2subroutine smoothsim(yt, ymiss, timevar, zt, ht,tt, rtv,qt,rqr, a1, p1, p1inf, &
3d, j, p, m, n, r,tol,rankp,ft,finf,kt,kinf,epshat,etahat,rt0,rt1,needeps)
4
5    implicit none
6
7    logical, intent(in) :: needeps
8    integer, intent(in) ::  p, m, n,r
9    integer, intent(inout) :: d, j,rankp
10    integer ::  t, i,tv
11    integer, intent(in), dimension(n,p) :: ymiss
12    integer, intent(in), dimension(5) :: timevar
13    double precision, intent(in), dimension(n,p) :: yt
14    double precision, intent(in), dimension(p,m,(n-1)*timevar(1)+1) :: zt
15    double precision, intent(in), dimension(p,p,(n-1)*timevar(2)+1) :: ht
16    double precision, intent(in), dimension(m,m,(n-1)*timevar(3)+1) :: tt
17    double precision, intent(in), dimension(m,r,(n-1)*timevar(4)+1) :: rtv
18    double precision, intent(in), dimension(r,r,(n-1)*timevar(5)+1) :: qt
19    double precision, dimension(m,m,(n-1)*max(timevar(4),timevar(5))+1) :: rqr
20    double precision, intent(in), dimension(m) :: a1
21    double precision, intent(in), dimension(m,m) ::  p1,p1inf
22    double precision, intent(inout), dimension(p,n) :: ft,finf
23    double precision, intent(inout), dimension(m,p,n) :: kt,kinf
24    double precision, dimension(p,n) :: vt
25    double precision, dimension(m) :: at
26    double precision, dimension(m,m) :: pt,pinf
27    double precision, intent(in) :: tol
28    double precision, dimension(m,m) :: im
29    double precision, intent(inout), dimension(r,n) :: etahat
30    double precision, intent(inout), dimension(p,n) :: epshat
31    double precision, intent(inout), dimension(m) :: rt0,rt1
32    double precision :: lik
33    double precision, external :: ddot
34
35    external dgemm, dsymm, dgemv, dsymv, dsyr, dsyr2, dger
36
37    tv = max(timevar(4),timevar(5))
38    lik = 0.0d0
39
40    j=0
41    d=0
42    pinf = p1inf
43    pt = p1
44    at = a1
45    if(rankp .GT. 0) then
46        diffuse: do while(d .LT. n .AND. rankp .GT. 0)
47            d = d+1
48            call dfilter1step(ymiss(d,:),yt(d,:),transpose(zt(:,:,(d-1)*timevar(1)+1)),ht(:,:,(d-1)*timevar(2)+1),&
49            tt(:,:,(d-1)*timevar(3)+1),rqr(:,:,(d-1)*tv+1),&
50            at,pt,vt(:,d),ft(:,d),kt(:,:,d),pinf,finf(:,d),kinf(:,:,d),rankp,lik,tol,0.0d0,p,m,j)
51        end do diffuse
52        if(rankp .EQ. 0 .AND. j .LT. p) then
53            call filter1step(ymiss(d,:),yt(d,:),transpose(zt(:,:,(d-1)*timevar(1)+1)),ht(:,:,(d-1)*timevar(2)+1),&
54            tt(:,:,(d-1)*timevar(3)+1),rqr(:,:,(d-1)*tv+1),&
55            at,pt,vt(:,d),ft(:,d),kt(:,:,d),lik,tol,0.0d0,p,m,j)
56        else
57            j = p
58        end if
59    end if
60
61    !Non-diffuse filtering continues from t=d+1, i=1
62
63    do t = d+1, n
64        call filter1step(ymiss(t,:),yt(t,:),transpose(zt(:,:,(t-1)*timevar(1)+1)),ht(:,:,(t-1)*timevar(2)+1),&
65        tt(:,:,(t-1)*timevar(3)+1),rqr(:,:,(t-1)*tv+1),&
66        at,pt,vt(:,t),ft(:,t),kt(:,:,t),lik,tol,0.0d0,p,m,0)
67    end do
68
69    !smoothing begins
70
71    im = 0.0d0
72    do i = 1, m
73        im(i,i) = 1.0d0
74    end do
75
76    rt0 = 0.0d0
77
78    do t = n, d+1, -1
79        call smooth1step(ymiss(t,:), transpose(zt(:,:,(t-1)*timevar(1)+1)), ht(:,:,(t-1)*timevar(2)+1), &
80        tt(:,:,(t-1)*timevar(3)+1), rtv(:,:,(t-1)*timevar(4)+1), qt(:,:,(t-1)*timevar(5)+1), vt(:,t), &
81        ft(:,t),kt(:,:,t), im,p,m,r,1,rt0,etahat(:,t),epshat(:,t),needeps)
82    end do
83
84    if(d .GT. 0) then
85        t = d
86        if(j .LT. p) then
87            call smooth1step(ymiss(t,:), transpose(zt(:,:,(t-1)*timevar(1)+1)), ht(:,:,(t-1)*timevar(2)+1), &
88            tt(:,:,(t-1)*timevar(3)+1), rtv(:,:,(t-1)*timevar(4)+1), qt(:,:,(t-1)*timevar(5)+1), vt(:,t), &
89            ft(:,t),kt(:,:,t), im,p,m,r,j+1,rt0,etahat(:,t),epshat(:,t),needeps)
90        end if
91        rt1 = 0.0d0
92        call dsmooth1step(ymiss(t,:), transpose(zt(:,:,(t-1)*timevar(1)+1)), ht(:,:,(t-1)*timevar(2)+1), &
93        tt(:,:,(t-1)*timevar(3)+1), rtv(:,:,(t-1)*timevar(4)+1), qt(:,:,(t-1)*timevar(5)+1), vt(:,t), &
94        ft(:,t),kt(:,:,t), im,p,m,r,j,rt0,rt1,finf(:,t),kinf(:,:,t),etahat(:,t),epshat(:,t),needeps)
95        do t=(d-1), 1, -1
96            call dsmooth1step(ymiss(t,:), transpose(zt(:,:,(t-1)*timevar(1)+1)), ht(:,:,(t-1)*timevar(2)+1), &
97            tt(:,:,(t-1)*timevar(3)+1), rtv(:,:,(t-1)*timevar(4)+1), qt(:,:,(t-1)*timevar(5)+1), vt(:,t), &
98            ft(:,t),kt(:,:,t), im,p,m,r,p,rt0,rt1,finf(:,t),kinf(:,:,t),etahat(:,t),epshat(:,t),needeps)
99        end do
100    end if
101
102end subroutine smoothsim
103