1!fast disturbance smoother for simulation
2subroutine smoothsimfast(yt, ymiss, timevar, zt, ht,tt, rtv,qt,a1, ft,kt,&
3finf, kinf, dt, jt, p, m, n,r,epshat,etahat,rt0,rt1,needeps)
4
5    implicit none
6
7    logical, intent(in) :: needeps
8    integer, intent(in) ::  p, m, r,n,dt,jt
9    integer ::  t, i
10    integer, intent(in), dimension(n,p) :: ymiss
11    integer, intent(in), dimension(5) :: timevar
12    double precision, intent(in), dimension(n,p) :: yt
13    double precision, intent(in), dimension(p,m,(n-1)*timevar(1)+1) :: zt
14    double precision, intent(in), dimension(p,p,(n-1)*timevar(2)+1) :: ht
15    double precision, intent(in), dimension(m,m,(n-1)*timevar(3)+1) :: tt
16    double precision, intent(in), dimension(m,r,(n-1)*timevar(4)+1) :: rtv
17    double precision, intent(in), dimension(r,r,(n-1)*timevar(5)+1) :: qt
18    double precision, intent(in), dimension(m) :: a1
19    double precision, intent(in), dimension(p,n) :: ft,finf
20    double precision, intent(in), dimension(m,p,n) :: kt,kinf
21    double precision, intent(inout), dimension(p,n) :: epshat
22    double precision, intent(inout), dimension(r,n) :: etahat
23    double precision, dimension(p,n) :: vt
24    double precision, dimension(m) :: at
25    double precision, dimension(m,m) :: im
26    double precision, intent(inout), dimension(m) :: rt0,rt1
27    double precision :: lik
28    lik = 0.0d0
29    at = a1
30    if(dt .GT. 0) then
31        !diffuse filtering begins
32        do t = 1, dt - 1
33            call dfilter1stepnv(ymiss(t,:),yt(t,:),transpose(zt(:,:,(t-1)*timevar(1)+1)),&
34            tt(:,:,(t-1)*timevar(3)+1),at,vt(:,t),ft(:,t),kt(:,:,t),&
35            finf(:,t),kinf(:,:,t),p,m,p,lik)
36        end do
37
38        t = dt
39        call dfilter1stepnv(ymiss(t,:),yt(t,:),transpose(zt(:,:,(t-1)*timevar(1)+1)),&
40        tt(:,:,(t-1)*timevar(3)+1),at,vt(:,t),ft(:,t),kt(:,:,t),&
41        finf(:,t),kinf(:,:,t),p,m,jt,lik)
42        !non-diffuse filtering begins
43        if(jt .LT. p) then
44            call filter1stepnv(ymiss(t,:),yt(t,:),transpose(zt(:,:,(t-1)*timevar(1)+1)),&
45            tt(:,:,(t-1)*timevar(3)+1),at,vt(:,t),ft(:,t),kt(:,:,t),p,m,jt,lik)
46        end if
47    end if
48    !Non-diffuse filtering continues from t=d+1, i=1
49    do t = dt + 1, n
50        call filter1stepnv(ymiss(t,:),yt(t,:),transpose(zt(:,:,(t-1)*timevar(1)+1)),&
51        tt(:,:,(t-1)*timevar(3)+1),at,vt(:,t),ft(:,t),kt(:,:,t),p,m,0,lik)
52    end do
53
54
55    !smoothing begins
56
57    im = 0.0d0
58    do i = 1, m
59        im(i,i) = 1.0d0
60    end do
61
62    rt0 = 0.0d0
63
64    do t = n, dt+1, -1
65        call smooth1step(ymiss(t,:), transpose(zt(:,:,(t-1)*timevar(1)+1)), ht(:,:,(t-1)*timevar(2)+1), &
66        tt(:,:,(t-1)*timevar(3)+1), rtv(:,:,(t-1)*timevar(4)+1), qt(:,:,(t-1)*timevar(5)+1), vt(:,t), &
67        ft(:,t),kt(:,:,t), im,p,m,r,1,rt0,etahat(:,t),epshat(:,t),needeps)
68    end do
69
70    if(dt .GT. 0) then
71        t = dt
72        if(jt .LT. p) then
73            call smooth1step(ymiss(t,:), transpose(zt(:,:,(t-1)*timevar(1)+1)), ht(:,:,(t-1)*timevar(2)+1), &
74            tt(:,:,(t-1)*timevar(3)+1), rtv(:,:,(t-1)*timevar(4)+1), qt(:,:,(t-1)*timevar(5)+1), vt(:,t), &
75            ft(:,t),kt(:,:,t), im,p,m,r,jt+1,rt0,etahat(:,t),epshat(:,t),needeps)
76        end if
77        rt1 = 0.0d0
78        call dsmooth1step(ymiss(t,:), transpose(zt(:,:,(t-1)*timevar(1)+1)), ht(:,:,(t-1)*timevar(2)+1), &
79        tt(:,:,(t-1)*timevar(3)+1), rtv(:,:,(t-1)*timevar(4)+1), qt(:,:,(t-1)*timevar(5)+1), vt(:,t), &
80        ft(:,t),kt(:,:,t), im,p,m,r,jt,rt0,rt1,finf(:,t),kinf(:,:,t),etahat(:,t),epshat(:,t),needeps)
81        do t = (dt - 1), 1, -1
82            call dsmooth1step(ymiss(t,:), transpose(zt(:,:,(t-1)*timevar(1)+1)), ht(:,:,(t-1)*timevar(2)+1), &
83            tt(:,:,(t-1)*timevar(3)+1), rtv(:,:,(t-1)*timevar(4)+1), qt(:,:,(t-1)*timevar(5)+1), vt(:,t), &
84            ft(:,t),kt(:,:,t), im,p,m,r,p,rt0,rt1,finf(:,t),kinf(:,:,t),etahat(:,t),epshat(:,t),needeps)
85        end do
86    end if
87
88end subroutine smoothsimfast
89