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