1!diffuse filtering for single time point 2subroutine dfilter1stepnv(ymiss, yt, zt, tt, at, vt, ft,kt,& 3finf,kinf,p,m,j,lik) 4 5 implicit none 6 7 integer, intent(in) :: p, m, j 8 integer :: i 9 integer, intent(in), dimension(p) :: ymiss 10 double precision, intent(in), dimension(p) :: yt 11 double precision, intent(in), dimension(m,p) :: zt 12 double precision, intent(in), dimension(m,m) :: tt 13 double precision, intent(inout), dimension(m) :: at 14 double precision, intent(inout), dimension(p) :: vt 15 double precision, intent(in), dimension(p) :: ft,finf 16 double precision, intent(in), dimension(m,p) :: kt,kinf 17 double precision, intent(inout) :: lik 18 double precision, dimension(m) :: ahelp 19 20 double precision, external :: ddot 21 external dgemv 22 23 do i = 1, j 24 if(ymiss(i).EQ.0) then 25 vt(i) = yt(i) - ddot(m,zt(:,i),1,at,1) 26 if (finf(i) .GT. 0.0d0) then 27 at = at + vt(i)/finf(i)*kinf(:,i) 28 lik = lik - 0.5d0*log(finf(i)) 29 else 30 if (ft(i) .GT. 0.0d0) then 31 at = at + vt(i)/ft(i)*kt(:,i) 32 lik = lik - 0.5d0*(log(ft(i)) + vt(i)**2/ft(i)) 33 end if 34 end if 35 end if 36 end do 37 if(j .EQ. p) then 38 call dgemv('n',m,m,1.0d0,tt,m,at,1,0.0d0,ahelp,1) 39 at = ahelp 40 end if 41end subroutine dfilter1stepnv 42 43 44!non-diffuse filtering for single time point 45subroutine filter1stepnv(ymiss, yt, zt, tt, at,vt, ft,kt,p,m,j,lik) 46 47 implicit none 48 49 integer, intent(in) :: p, m,j 50 integer :: i 51 integer, intent(in), dimension(p) :: ymiss 52 double precision, intent(in), dimension(p) :: yt 53 double precision, intent(in), dimension(m,p) :: zt 54 double precision, intent(in), dimension(m,m) :: tt 55 double precision, intent(inout), dimension(m) :: at 56 double precision, intent(in), dimension(p) :: ft 57 double precision, intent(inout), dimension(p) :: vt 58 double precision, intent(in), dimension(m,p) :: kt 59 double precision, intent(inout) :: lik 60 double precision, dimension(m) :: ahelp 61 62 double precision, external :: ddot 63 64 external dgemv 65 66 do i = j+1, p 67 if(ymiss(i).EQ.0) then 68 vt(i) = yt(i) - ddot(m,zt(:,i),1,at,1) 69 if (ft(i) .GT. 0.0d0) then 70 at = at + vt(i)/ft(i)*kt(:,i) 71 lik = lik - 0.5d0*(log(ft(i)) + vt(i)**2/ft(i)) 72 end if 73 end if 74 end do 75 76 call dgemv('n',m,m,1.0d0,tt,m,at,1,0.0d0,ahelp,1) 77 at = ahelp 78 79end subroutine filter1stepnv 80