    PROGRAM g13bjfe

!      G13BJF Example Program Text

!      Mark 23 Release. NAG Copyright 2011.

!      .. Use Statements ..
       USE nag_library, ONLY : g13bjf, nag_wp, x04caf
!      .. Implicit None Statement ..
       IMPLICIT NONE
!      .. Parameters ..
       INTEGER, PARAMETER              :: nin = 5, nout = 6
!      .. Local Scalars ..
       INTEGER                         :: dp, i, ifail, imwa, isttf, iwa, kfc, &
                                          kzef, ldparx, ldxxy, mx, n, ncf,     &
                                          ncg, nch, nci, nev, nfv, nis, npara, &
                                          nparx, nser, nsttf, qp, qx, smx
!      .. Local Arrays ..
       REAL (KIND=nag_wp), ALLOCATABLE :: fsd(:), fva(:), para(:), parx(:,:),  &
                                          rmsxy(:), sttf(:), wa(:), xxy(:,:)
       INTEGER                         :: mr(7)
       INTEGER, ALLOCATABLE            :: mrx(:,:), mt(:,:), mwa(:)
!      .. Intrinsic Functions ..
       INTRINSIC                          max, sum
!      .. Executable Statements ..
       WRITE (nout,*) 'G13BJF Example Program Results'
       WRITE (nout,*)

!      Skip heading in data file
       READ (nin,*)

!      Read in problem size
       READ (nin,*) kfc, nev, nfv, nser, kzef

!      Number of input series
       nis = nser - 1

       ALLOCATE (mt(4,nser))

!      Read in the orders for the output noise
       READ (nin,*) mr(1:7)

!      Read in transfer function
       DO i = 1, nis
          READ (nin,*) mt(1:4,i)
       END DO

!      Calculate NPARA
       npara = 0
       DO i = 1, nis
          npara = npara + mt(2,i) + mt(3,i)
       END DO
       npara = npara + mr(1) + mr(3) + mr(4) + mr(6) + nser

!      Calculate array sizes
       n = nev + nfv
       ldxxy = n
       ncf = 0
       DO i = 1, nis
          IF (mt(4,i)>1) THEN
             ncf = sum(mt(1:3,i))
          END IF
       END DO
       isttf = mr(4)*mr(7) + mr(2) + mr(5)*mr(7) + mr(3) + &
          max(mr(1),mr(6)*mr(7)) + ncf
       qp = mr(3) + mr(6)*mr(7)
       dp = mr(2) + mr(5)*mr(7)
       smx = 0
       qx = qp
       nci = nser
       DO i = 1, nis
          IF (mt(4,i)==3) THEN
             mx = max(mt(1,i)+mt(2,i),mt(3,i))
             nci = nci + 1
          ELSE
             mx = 0
          END IF
          IF (mt(3,i)>0) THEN
             nci = nci + 1
          END IF
          smx = smx + mx
          qx = max(qx,mx)
       END DO
       ncg = npara + qx + smx
       nch = dp + 6*qx + nev
       IF (qx>0) THEN
          nci = nci + 1
       END IF
       IF (mr(1)>0) THEN
          nci = nci + 1
       END IF
       IF (mr(3)>0) THEN
          nci = nci + 1
       END IF
       IF (mr(4)>0) THEN
          nci = nci + 1
       END IF
       IF (mr(6)>0) THEN
          nci = nci + 1
       END IF
       iwa = 2*(ncg**2) + nch*(nci+4)
       imwa = 16*nser + 7*ncg + 3*npara + 27
       ALLOCATE (para(npara),xxy(ldxxy,nser),rmsxy(nser),mrx(7,nser),fva(nfv), &
          fsd(nfv),sttf(isttf),wa(iwa),mwa(imwa))

!      Read in multi-input model parameters
       READ (nin,*) para(1:npara)

!      Read in the observed values for the input and output series
       READ (nin,*) (xxy(i,1:nser),i=1,nev)

!      Read in the future values for the input series
       READ (nin,*) (xxy(nev+i,1:nis),i=1,nfv)

       IF (nis>=1) THEN
!         Read in residual variance of input series
          READ (nin,*) rmsxy(1:nis)

!         Read in orders for input series ARIMA where available
!         (i.e. where residual variance is not zero)
          ldparx = 0
          DO i = 1, nis
             IF (rmsxy(i)/=0.0E0_nag_wp) THEN
                READ (nin,*) mrx(1:7,i)
                nparx = mrx(1,i) + mrx(3,i) + mrx(4,i) + mrx(6,i)
                ldparx = max(ldparx,nparx)
             ELSE
                mrx(1:7,i) = 0
             END IF
          END DO
       ELSE
!         No input series
          ldparx = 1
       END IF

       ALLOCATE (parx(ldparx,nser))

!      Read in parameters for each input series ARIMA         
       IF (nis>0) THEN
          DO i = 1, nis
             IF (rmsxy(i)/=0.0E0_nag_wp) THEN
                nparx = mrx(1,i) + mrx(3,i) + mrx(4,i) + mrx(6,i)
                IF (nparx>0) THEN
                   READ (nin,*) parx(1:nparx,i)
                END IF
             END IF
          END DO
       END IF

       ifail = 0
       CALL g13bjf(mr,nser,mt,para,npara,kfc,nev,nfv,xxy,ldxxy,kzef,rmsxy,mrx, &
          parx,ldparx,fva,fsd,sttf,isttf,nsttf,wa,iwa,mwa,imwa,ifail)

!      Display results
       WRITE (nout,99999) 'After processing', nev, ' sets of observations'
       WRITE (nout,99998) nsttf, ' values of the state set are derived'
       WRITE (nout,*)
       WRITE (nout,99997) sttf(1:nsttf)
       WRITE (nout,*)
       WRITE (nout,*) 'The residual mean square for the output'
       WRITE (nout,99996) 'series is also derived and its value is', &
          rmsxy(nser)
       WRITE (nout,*)
       WRITE (nout,*) 'The forecast values and their standard errors are'
       WRITE (nout,*)
       WRITE (nout,*) '   I       FVA       FSD'
       WRITE (nout,*)
       WRITE (nout,99995) (i,fva(i),fsd(i),i=1,nfv)
       WRITE (nout,*)
       FLUSH (nout)
       ifail = 0
       CALL x04caf('General',' ',n,nser,xxy,ldxxy, &
          'The values of z(t) and n(t) are',ifail)
       WRITE (nout,99994) 'The first ', nis, &
          ' columns hold the z(t) and the last column the n(t)'

99999  FORMAT (1X,A,I3,A)
99998  FORMAT (1X,I3,A)
99997  FORMAT (1X,6F10.4)
99996  FORMAT (1X,A,F10.4)
99995  FORMAT (1X,I4,F10.3,F10.4)
99994  FORMAT (1X,A,I0,A)
    END PROGRAM g13bjfe
