diff --git a/prep_syndata/calc_error_cov.f90 b/prep_syndata/calc_error_cov.f90
new file mode 100644
index 0000000000000000000000000000000000000000..d091d2792caea26d62bb703e95dbf2b6ecd53e52
--- /dev/null
+++ b/prep_syndata/calc_error_cov.f90
@@ -0,0 +1,209 @@
+!---------------------------------------------------------------------------------------
+! FLEXINVERT: calc_error_cov
+!---------------------------------------------------------------------------------------
+! FLEXINVERT is free software: you can redistribute it and/or modify
+! it under the terms of the GNU General Public License as published by
+! the Free Software Foundation, either version 3 of the License, or
+! (at your option) any later version.
+!
+! FLEXINVERT is distributed in the hope that it will be useful,
+! but WITHOUT ANY WARRANTY; without even the implied warranty of
+! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+! GNU General Public License for more details.
+!
+! You should have received a copy of the GNU General Public License
+! along with FLEXINVERT. If not, see <http://www.gnu.org/licenses/>.
+!
+! Copyright 2017, Rona Thompson
+!---------------------------------------------------------------------------------------
+!
+!> calc_error_cov
+!!
+!! Purpose: Calculates the prior error covariance matrix for the state variables.
+!!
+!! Interface:
+!!
+!! Inputs
+!! config - configuration settings data structure
+!! files - file settings data structure
+!!
+!! Externals
+!! dsyevr (LAPACK routine)
+!!
+!---------------------------------------------------------------------------------------
+
+subroutine calc_error_cov(config, files, state)
+
+ use mod_var
+ use mod_settings
+ use mod_tools
+
+ implicit none
+
+ type (config_t), intent (in) :: config
+ type (files_t), intent (in) :: files
+ real, dimension(nvar), intent (in out) :: state
+
+ character(len=max_path_len) :: filename
+ character(len=max_path_len) :: rowfmt
+ character(len=max_name_len) :: varname, lonname, latname, timename
+ integer :: i, j, k, it, jt, nt, num, ix, jy, nb, nread
+ real(kind=8) :: jd
+ real, dimension(ndvar) :: xerr
+ real(kind=8), dimension(ndvar,ndvar) :: cov
+ real, dimension(nbox,nbox) :: corr
+ real, dimension(:,:,:), allocatable :: tmp
+ integer :: lwork, liwork
+ integer, dimension(:), allocatable :: isuppz
+ integer, dimension(:), allocatable :: iwork
+ real(kind=8), dimension(:), allocatable :: work
+ real(kind=8), dimension(:,:), allocatable :: zwork
+ real(kind=8), dimension(:), allocatable :: evaltmp
+ real(kind=8), dimension(:,:), allocatable :: evectmp
+ integer :: ierr, info, neval, ind
+ real, parameter :: smallnum=1.e-15
+ real :: area
+ logical :: lexist
+
+ ! read error from analysis
+ ! ------------------------
+
+ inquire(file=trim(files%path_output)//'analysis.nc',exist=lexist)
+ if ( lexist ) then
+ filename = trim(files%path_output)//'analysis.nc'
+ lonname = 'longitude'
+ latname = 'latitude'
+ timename = 'time'
+ jd = juldatei
+ ! prior error
+ varname = 'epri'
+ nread = 1
+ allocate( tmp(nxregrid,nyregrid,1) )
+ call read_ncdf(filename, &
+ varname, &
+ lonname, &
+ latname, &
+ timename,&
+ rllx, rlly, &
+ nxregrid, nyregrid, &
+ jd, nread, num, &
+ tmp)
+ tmp = tmp*numscale
+ xerr(:) = 0.
+ do jy = 1, nyregrid
+ area = grid_area(reg_lat(jy)+0.5*rdy, rdy, rdx)
+ do ix = 1, nxregrid
+ xerr(nbox_xy(ix,jy)) = xerr(nbox_xy(ix,jy)) + tmp(ix,jy,1)*area
+ end do
+ end do
+ xerr(:) = xerr(:)/area_box(:)
+ deallocate(tmp)
+ ! prior flux
+ varname = 'fpri'
+ nread = ntstate
+ allocate( tmp(nxregrid,nyregrid,ntstate) )
+ call read_ncdf(filename, &
+ varname, &
+ lonname, &
+ latname, &
+ timename,&
+ rllx, rlly, &
+ nxregrid, nyregrid, &
+ jd, nread, num, &
+ tmp)
+ tmp = tmp*numscale
+ do nt = 1, ntstate
+ do jy = 1, nyregrid
+ area = grid_area(reg_lat(jy)+0.5*rdy, rdy, rdx)
+ do ix = 1, nxregrid
+ state((nt-1)*nbox+nbox_xy(ix,jy)) = state((nt-1)*nbox+nbox_xy(ix,jy)) + tmp(ix,jy,nt)*area
+ end do
+ end do
+ state((nt-1)*nbox+1:nt*nbox) = state((nt-1)*nbox+1:nt*nbox)/area_box(:)
+ end do
+ deallocate(tmp)
+ else
+ write(logid,*) 'ERROR: file analysis.nc not found'
+ stop
+ endif
+
+ ! read correlations
+ ! -----------------
+
+ inquire(file=trim(files%path_output)//'correl.txt',exist=lexist)
+ if ( lexist ) then
+ write(logid,*) 'Reading spatial correlation from correl.txt'
+ open(100,file=trim(files%path_output)//'correl.txt',action='read',status='old',iostat=ierr)
+ write(rowfmt,'(A,I6,A)') '(',nbox,'(E11.4,1X))'
+ print*, rowfmt
+ do nb = 1, nbox
+ read(100,fmt=rowfmt) corr(nb,:)
+ end do
+ close(100)
+ else
+ write(logid,*) 'ERROR: cannot find correl.txt'
+ stop
+ endif
+
+ ! calculate error covariance
+ ! --------------------------
+
+ ! spatial error covariance one time step
+ cov(:,:) = 0d0
+ ! sub-daily resolution (no correlation between sub-daily time steps)
+ do nt = 1, 1
+ do j = 1, nbox
+ do i = 1, nbox
+ cov((nt-1)*nbox+i,(nt-1)*nbox+j) = dble(corr(i,j))*dble(xerr((nt-1)*nbox+i))*dble(xerr((nt-1)*nbox+j))
+ end do
+ end do
+ end do
+
+ ! calculate eigendecomposition of cov
+ ! -----------------------------------
+
+ ! B = U*D*U^(T)
+ lwork = 26*ndvar
+ liwork = 10*ndvar
+ allocate( work(lwork), stat=ierr )
+ if ( ierr.ne.0 ) stop 'ERROR error_cov: not enough memory'
+ allocate( iwork(liwork), stat=ierr )
+ if ( ierr.ne.0 ) stop 'ERROR error_cov: not enough memory'
+ allocate( isuppz(2*ndvar), stat=ierr )
+ if ( ierr.ne.0 ) stop 'ERROR error_cov: not enough memory'
+ allocate( evectmp(ndvar,ndvar), stat=ierr )
+ if ( ierr.ne.0 ) stop 'ERROR error_cov: not enough memory'
+ allocate( evaltmp(ndvar), stat=ierr )
+ if ( ierr.ne.0 ) stop 'ERROR error_cov: not enough memory'
+ allocate( zwork(ndvar,ndvar), stat=ierr )
+ if ( ierr.ne.0 ) stop 'ERROR error_cov: not enough memory'
+
+ zwork = cov
+ call dsyevr('V','A','L',ndvar,zwork,ndvar,1.,1.,1,1,0.,neval,evaltmp,evectmp,ndvar,isuppz,work,lwork,iwork,liwork,info)
+ if ( info.ne.0 ) then
+ write(logid,*) 'ERROR error_cov: eigenvalue decomposition'
+ stop
+ endif
+
+ ! truncate eigenvalues
+ ! eigenvalues are sorted in ascending order
+ ind = minloc(evaltmp, dim=1, mask=evaltmp.ge.(trunc*maxval(evaltmp)))
+ write(logid,*) 'Truncating eigenvalues at: ',trunc*maxval(evaltmp)
+ neig = ndvar - ind + 1
+ write(logid,*) 'Number of retained eigenvalues: ',neig
+
+ allocate( evals(neig) )
+ allocate( evecs(ndvar, neig) )
+ evals = evaltmp(ind:ndvar)
+ evecs = evectmp(:,ind:ndvar)
+ deallocate( evaltmp )
+ deallocate( evectmp)
+
+ deallocate(work)
+ deallocate(zwork)
+ deallocate(iwork)
+ deallocate(isuppz)
+
+
+end subroutine calc_error_cov
+