added calc_error_cov for prep_synthdata

prep_syndata/calc_error_cov.f90

+!---------------------------------------------------------------------------------------
+! 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
+