From 52643eac60442d2cdef9406d4ab6161ecb5c6fc7 Mon Sep 17 00:00:00 2001
From: ronesy <rlt@nilu.no>
Date: Mon, 3 Mar 2025 13:57:12 +0100
Subject: [PATCH] Changes to include temporal correlation between scalars of
initial mixing ratios
---
source/congrad.f90 | 20 ++++-
source/error_cov.f90 | 139 +++++++++++++++++-------------
source/initialize.f90 | 3 +
source/mod_analytic.f90 | 73 +++++++++++++---
source/mod_covar.f90 | 171 ++++++++++++++++++++++++++++++++++---
source/mod_perturb.f90 | 15 +++-
source/mod_save.f90 | 178 ++++++++++++++++++++++++++++-----------
source/mod_settings.f90 | 45 ++++++----
source/mod_transform.f90 | 119 ++++++++++++++++++++------
source/simulate.f90 | 3 +
10 files changed, 589 insertions(+), 177 deletions(-)
diff --git a/source/congrad.f90 b/source/congrad.f90
index 8c1373f..4ef674b 100644
--- a/source/congrad.f90
+++ b/source/congrad.f90
@@ -93,6 +93,7 @@ subroutine congrad(iter, grad, files, config, fluxes, obs, states, covar)
real(kind=8), dimension(maxiter+1,1) :: zwork
real(kind=8), dimension(neig) :: tmp0
real(kind=8), dimension(ndvar) :: tmp1
+ real(kind=8), dimension(ntcini*ncini) :: tmp2
real(kind=8), dimension(nvar) :: tmp
real(kind=8) :: dla
real(kind=8) :: zgnorm ! initial gradient norm
@@ -385,16 +386,19 @@ subroutine congrad(iter, grad, files, config, fluxes, obs, states, covar)
close(100)
! calculate inverse of Hessian
+ ! xerr^2 = xerr0^2 - sum[(1-1/eval[j])(B^(1/2)*evec(j,))^2]
+ ! calculate sum[(1-1/eval[j])(B^(1/2)*evec(j,))^2]
+ ! first calculate pevecs = (1-1/eval[j])^(1/2)*evec(j,)
allocate( pevecs(ngood,nvar) )
pevecs(:,:) = 0d0
do jk = 1, ngood
pevecs(jk,:) = pvcglev(:,jk)*sqrt(1d0-1d0/ptheta(jk))
enddo
- ! xerr^2 = xerr0^2 - sum[(1-1/eval[j])(B^(1/2)*evec(j,))^2]
- ! calculate sum[(1-1/eval[j])(B^(1/2)*evec(j,))^2]
+ ! second calculate B^(1/2)*pevecs
tmp(:) = 0d0
do jk = 1, ngood
+ ! fluxes
do jt = 1, ntstate
do k = 1, ndvar
do j = 1, neig
@@ -409,8 +413,18 @@ subroutine congrad(iter, grad, files, config, fluxes, obs, states, covar)
end do
end do
end do
+ ! scalars of initial mixing ratios
if ( config%opt_cini ) then
- tmp(npvar+1:nvar) = tmp(npvar+1:nvar) + (config%cinierr * pevecs(jk,npvar+1:nvar))**2
+ do jt = 1, ntcini
+ do k = 1, ncini
+ do j = 1, ntcini
+ tmp2((j-1)*ncini+1:j*ncini) = covar%sqcortcini(jt,j)*covar%sqcorcini(k,:)
+ end do
+ tmp(npvar+(jt-1)*ncini+k) = tmp(npvar+(jt-1)*ncini+k) + &
+ dot_product(tmp2,pevecs(jk,npvar+1:nvar))**2
+ end do
+ end do
+! tmp(npvar+1:nvar) = tmp(npvar+1:nvar) + (config%cinierr * pevecs(jk,npvar+1:nvar))**2
endif
end do
! calculate last part
diff --git a/source/error_cov.f90 b/source/error_cov.f90
index 2ab3dc6..a22856e 100644
--- a/source/error_cov.f90
+++ b/source/error_cov.f90
@@ -62,8 +62,8 @@ subroutine error_cov(config, files, states, covar, corr, xerr)
real(kind=8), dimension(ndvar,ndvar) :: cov
real(kind=8), dimension(ntstate) :: tmp
real(kind=8), dimension(nvar) :: var
- real, dimension(ntstate) :: tsteps
- real, dimension(ntstate,ntstate) :: tmat, dtmat
+ real, dimension(:), allocatable :: tsteps
+ real, dimension(:,:), allocatable :: tmat, dtmat
real(kind=8), dimension(ntstate,ntstate) :: cort, icort, sqcort, isqcort
integer :: lwork, liwork
integer, dimension(:), allocatable :: isuppz
@@ -75,6 +75,9 @@ subroutine error_cov(config, files, states, covar, corr, xerr)
integer :: ierr, info, neval, ind
real, parameter :: smallnum=1.e-15
integer :: timei, timef, countrate
+ real(kind=8), dimension(:,:), allocatable :: imat, sqmat, isqmat
+ real(kind=8), dimension(ntcini,ntcini) :: cortcini
+ real(kind=8), dimension(ncini,ncini) :: corcini
! calculate error covariance
! --------------------------
@@ -207,6 +210,9 @@ subroutine error_cov(config, files, states, covar, corr, xerr)
! calculate temporal correlation
! ------------------------------
+ allocate( tsteps(ntstate) )
+ allocate( tmat(ntstate,ntstate) )
+ allocate( dtmat(ntstate,ntstate) )
do i = 1, ntstate
tsteps(i) = i * statres
end do
@@ -220,71 +226,72 @@ subroutine error_cov(config, files, states, covar, corr, xerr)
cort(j,i) = exp(-1d0*dble(dtmat(j,i)/config%sigmatime))
end do
end do
+ deallocate( tsteps, tmat, dtmat )
- ! calculate eigendecomposition of cort
- ! ------------------------------------
+ ! calculate matrix inverse, square root and inverse square root
+ allocate( imat(ntstate,ntstate), sqmat(ntstate,ntstate), isqmat(ntstate,ntstate) )
+ call matrix_oper(ntstate, cort, imat, sqmat, isqmat)
+ covar%cort = cort
+ covar%icort = imat
+ covar%sqcort = sqmat
+ covar%isqcort = isqmat
+ deallocate( imat, sqmat, isqmat )
- ! cort = U*D*U^(T)
- lwork = 26*ntstate
- liwork = 10*ntstate
- 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*ntstate), stat=ierr )
- if ( ierr.ne.0 ) stop 'ERROR error_cov: not enough memory'
- allocate( evecs(ntstate,ntstate), stat=ierr )
- if ( ierr.ne.0 ) stop 'ERROR error_cov: not enough memory'
- allocate( evals(ntstate), stat=ierr )
- if ( ierr.ne.0 ) stop 'ERROR error_cov: not enough memory'
- allocate( zwork(ntstate,ntstate), stat=ierr )
- if ( ierr.ne.0 ) stop 'ERROR error_cov: not enough memory'
+ ! scalars of initial mixing ratios
+ ! --------------------------------
- zwork = cort
- call dsyevr('V','A','L',ntstate,zwork,ntstate,1.,1.,1,1,0.,neval,evals,evecs,ntstate,isuppz,work,lwork,iwork,liwork,info)
- if ( info.ne.0 ) then
- write(logid,*) 'ERROR error_cov: eigenvalue decomposition'
- stop
- endif
+ if ( config%opt_cini ) then
- ! inverse of cort
- icort(:,:) = 0d0
- do i = 1, ntstate
- do j = 1, ntstate
- tmp(j) = evecs(i,j)*1d0/evals(j)
+ allocate( tsteps(ntcini) )
+ allocate( tmat(ntcini,ntcini) )
+ allocate( dtmat(ntcini,ntcini) )
+ do i = 1, ntcini
+ tsteps(i) = i * cinires
end do
- icort(i,:) = matmul(tmp,transpose(evecs))
- end do
-
- ! square root of cort
- sqcort(:,:) = 0d0
- do i = 1, ntstate
- do j = 1, ntstate
- tmp(j) = evecs(i,j)*sqrt(evals(j))
+ tmat = matmul(reshape(tsteps,(/ntcini,1/)),reshape((tsteps/tsteps),(/1,ntcini/)))
+ dtmat = abs(tmat - transpose(tmat))
+ print*, 'error_cov: dtmat = ',dtmat
+
+ ! temporal correlation matrix
+ cortcini(:,:) = 0d0
+ do i = 1, ntcini
+ do j = i, ntcini
+ cortcini(i,j) = exp(-1d0*dble(dtmat(i,j)/config%sigmatime_cini))
+ cortcini(j,i) = exp(-1d0*dble(dtmat(j,i)/config%sigmatime_cini))
+ end do
end do
- sqcort(i,:) = matmul(tmp,transpose(evecs))
- end do
-
- ! inverse of square root of cort
- isqcort(:,:) = 0d0
- do i = 1, ntstate
- do j = 1, ntstate
- tmp(j) = evecs(i,j)*1d0/sqrt(evals(j))
+ deallocate( tsteps, tmat, dtmat )
+
+ ! calculate matrix inverse, square root and inverse square root
+ allocate( imat(ntcini,ntcini), sqmat(ntcini,ntcini), isqmat(ntcini,ntcini) )
+ call matrix_oper(ntcini, cortcini, imat, sqmat, isqmat)
+ covar%cortcini = cortcini
+ covar%icortcini = imat
+ covar%sqcortcini = sqmat
+ covar%isqcortcini = isqmat
+ deallocate( imat, sqmat, isqmat )
+
+ ! covariance matrix for one time step
+ ! currently without spatial correlation
+ corcini(:,:) = 0d0
+ do i = 1, ncini
+ if ( config%lognormal ) then
+ corcini(i,i) = (log(1.+config%cinierr))**2
+ else
+ corcini(i,i) = config%cinierr**2
+ endif
end do
- isqcort(i,:) = matmul(tmp,transpose(evecs))
- end do
- covar%cort = cort
- covar%icort = icort
- covar%sqcort = sqcort
- covar%isqcort = isqcort
+ ! calculate matrix inverse, square root and inverse square root
+ allocate( imat(ncini,ncini), sqmat(ncini,ncini), isqmat(ncini,ncini) )
+ call matrix_oper(ncini, corcini, imat, sqmat, isqmat)
+ covar%corcini = corcini
+ covar%icorcini = imat
+ covar%sqcorcini = sqmat
+ covar%isqcorcini = isqmat
+ deallocate( imat, sqmat, isqmat )
- deallocate(work)
- deallocate(zwork)
- deallocate(iwork)
- deallocate(isuppz)
- deallocate(evals)
- deallocate(evecs)
+ endif
! prior error
! -----------
@@ -347,6 +354,22 @@ subroutine error_cov(config, files, states, covar, corr, xerr)
write(100,fmt=rowfmt) covar%sqcort(i,:)
end do
close(100)
+ if ( config%opt_cini ) then
+ filename = trim(files%path_output)//'cortcini.txt'
+ write(rowfmt,'(A,I6,A)') '(',ntcini,'(E11.4,1X))'
+ open(100,file=trim(filename),status='replace',action='write',iostat=ierr)
+ do i = 1, ntcini
+ write(100,fmt=rowfmt) covar%cortcini(i,:)
+ end do
+ close(100)
+ filename = trim(files%path_output)//'corcini.txt'
+ write(rowfmt,'(A,I6,A)') '(',ncini,'(E11.4,1X))'
+ open(100,file=trim(filename),status='replace',action='write',iostat=ierr)
+ do i = 1, ncini
+ write(100,fmt=rowfmt) covar%corcini(i,:)
+ end do
+ close(100)
+ endif
! verbose output
! --------------
diff --git a/source/initialize.f90 b/source/initialize.f90
index f0b7b1a..bab9330 100644
--- a/source/initialize.f90
+++ b/source/initialize.f90
@@ -386,6 +386,9 @@ subroutine initialize(files, config)
write(logid,*) 'nday: ',nday
write(logid,*) 'statres: ',statres
write(logid,*) 'ntstate: ',ntstate
+ if ( config%opt_cini ) then
+ write(logid,*) 'ntcini: ',ntcini
+ endif
if ( trim(config%spec).eq.'co2' ) then
write(logid,*) 'nd_nee: ',nd_nee
write(logid,*) 'nt_nee: ',nt_nee
diff --git a/source/mod_analytic.f90 b/source/mod_analytic.f90
index 146100e..c838dbc 100644
--- a/source/mod_analytic.f90
+++ b/source/mod_analytic.f90
@@ -153,6 +153,7 @@ module mod_analytic
stop
endif
zwork(:,:) = 0d0
+ ! fluxes
do jt = 1, ntstate
do k = 1, ndvar
do j = 1, neig
@@ -169,12 +170,22 @@ module mod_analytic
end do
end do
end do
+ ! scalars of initial mixing ratios
if ( config%opt_cini ) then
zwork(:,npvar+1:nvar) = matmul(nwork(:,:),dble(hmat(:,npvar+1:nvar)))
zwork(npvar+1:nvar,:) = matmul(nwork(npvar+1:nvar,:),dble(hmat(:,:)))
- do it = 1, ntcini*ncini
- zwork(npvar+it,npvar+it) = zwork(npvar+it,npvar+it) + 1./config%cinierr**2
+ do it = 1, ntcini
+ do k = 1, ncini
+ do jt = 1, ntcini
+ zwork(npvar+(it-1)*ncini+k,npvar+(jt-1)*ncini+1:npvar+jt*ncini) = &
+ zwork(npvar+(it-1)*ncini+k,npvar+(jt-1)*ncini+1:npvar+jt*ncini) + &
+ covar%icortcini(it,jt)*covar%icorcini(k,:)
+ end do
+ end do
end do
+! do it = 1, ntcini*ncini
+! zwork(npvar+it,npvar+it) = zwork(npvar+it,npvar+it) + 1./config%cinierr**2
+! end do
endif
! verbose only
@@ -192,19 +203,30 @@ module mod_analytic
stop
endif
! H^TR^(-1)H + B^(-1)
+ ! fluxes
do jt = 1, ntstate
do it = 1, ntstate
zwork_dir((jt-1)*ndvar+1:jt*ndvar,(it-1)*ndvar+1:it*ndvar) = &
matmul(nwork((jt-1)*ndvar+1:jt*ndvar,:),dble(hmat(:,(it-1)*ndvar+1:it*ndvar))) + covar%icort(jt,it)*icovb
end do
end do
+ ! scalars of initial mixing ratios
if ( config%opt_cini ) then
zwork_dir(:,npvar+1:nvar) = matmul(nwork(:,:),dble(hmat(:,npvar+1:nvar)))
zwork_dir(npvar+1:nvar,:) = matmul(nwork(npvar+1:nvar,:),dble(hmat(:,:)))
- do it = 1, ntcini*ncini
- zwork_dir(npvar+it,npvar+it) = zwork_dir(npvar+it,npvar+it) + 1./config%cinierr**2
+ do it = 1, ntcini
+ do k = 1, ncini
+ do jt = 1, ntcini
+ zwork_dir(npvar+(it-1)*ncini+k,npvar+(jt-1)*ncini+1:npvar+jt*ncini) = &
+ zwork_dir(npvar+(it-1)*ncini+k,npvar+(jt-1)*ncini+1:npvar+jt*ncini) + &
+ covar%icortcini(it,jt)*covar%icorcini(k,:)
+ end do
+ end do
end do
- endif
+! do it = 1, ntcini*ncini
+! zwork_dir(npvar+it,npvar+it) = zwork_dir(npvar+it,npvar+it) + 1./config%cinierr**2
+! end do
+ endif
filename = trim(files%path_output)//'zwork.txt'
open(100,file=trim(filename),status='replace',action='write',iostat=ierr)
write(rowfmt,'(A,I6,A)') '(',nvar,'(E11.4,1X))'
@@ -354,6 +376,8 @@ module mod_analytic
real(kind=8), dimension(:,:), allocatable :: imwork
real(kind=8), dimension(neig) :: tmp
real(kind=8), dimension(1,ndvar) :: tmp1
+ real(kind=8), dimension(1,ntcini*ncini) :: tmp2
+ real(kind=8), dimension(ntcini*ncini,ntcini*ncini) :: tmp3
integer :: i, j, k, it, jt
integer :: info, ierr
@@ -371,6 +395,7 @@ module mod_analytic
stop
endif
zwork(:,:) = 0d0
+ ! fluxes
do jt = 1, ntstate
do k = 1, ndvar
do j = 1, neig
@@ -387,8 +412,18 @@ module mod_analytic
end do
end do
end do
+ ! scalars of initial mixing ratio
if ( config%opt_cini ) then
- zwork(npvar+1:nvar,:) = config%cinierr**2 * transpose(dble(hmat(:,npvar+1:nvar)))
+ do it = 1, ntcini
+ do k = 1, ncini
+ do jt = 1, ntcini
+ ! row vector of B for scalars
+ tmp2(1,(jt-1)*ncini+1:jt*ncini) = covar%cortcini(it,jt)*covar%corcini(k,:)
+ end do
+ zwork(npvar+(it-1)*ncini+k,:) = matmul(tmp2(1,:),transpose(dble(hmat(:,npvar+1:nvar))))
+ end do
+ end do
+! zwork(npvar+1:nvar,:) = config%cinierr**2 * transpose(dble(hmat(:,npvar+1:nvar)))
endif
! verbose only
@@ -401,6 +436,7 @@ module mod_analytic
end do
covb = matmul(tmp0,transpose(covar%evecs))
! BH^(T)
+ ! fluxes
allocate ( zwork_dir(nvar,nobs), stat=ierr )
if ( ierr.ne.0 ) then
write(logid,*) 'ERROR invert_nobs: not enough memory'
@@ -413,8 +449,15 @@ module mod_analytic
matmul(covar%cort(jt,it)*covb,transpose(dble(hmat(:,(it-1)*ndvar+1:it*ndvar))))
end do
end do
+ ! scalars of initial mixing ratios
if ( config%opt_cini ) then
- zwork_dir(npvar+1:nvar,:) = config%cinierr**2 * transpose(dble(hmat(:,npvar+1:nvar)))
+ do it = 1, ntcini
+ do jt = 1, ntcini
+ tmp3((it-1)*ncini+1:it*ncini,(jt-1)*ncini+1:jt*ncini) = covar%cortcini(it,jt)*covar%corcini(:,:)
+ end do
+ end do
+ zwork_dir(npvar+1:nvar,:) = matmul(tmp3(:,:),transpose(dble(hmat(:,npvar+1:nvar))))
+! zwork_dir(npvar+1:nvar,:) = config%cinierr**2 * transpose(dble(hmat(:,npvar+1:nvar)))
endif
filename = trim(files%path_output)//'zwork.txt'
open(100,file=trim(filename),status='replace',action='write',iostat=ierr)
@@ -500,6 +543,7 @@ module mod_analytic
stop
endif
cova(:,:) = 0d0
+ ! fluxes
do jt = 1, ntstate
do k = 1, ndvar
do j = 1, neig
@@ -517,11 +561,20 @@ module mod_analytic
end do
end do
if ( config%opt_cini ) then
- cova(1:npvar,npvar+1:nvar) = -1d0*matmul(gain(1:npvar,:),transpose(zwork(npvar+1:nvar,:)))
+ ! subtract remaining part of GHB
+ cova(1:npvar,npvar+1:nvar) = -1d0*matmul(gain(1:npvar,:),transpose(zwork(npvar+1:nvar,:)))
cova(npvar+1:nvar,:) = -1d0*matmul(gain(npvar+1:nvar,:),transpose(zwork(:,:)))
- do i = 1, ntcini*ncini
- cova(npvar+i,npvar+i) = config%cinierr**2 + cova(npvar+i,npvar+i)
+ ! add B for scalars of initial mixing ratios
+ do it = 1, ntcini
+ do jt = 1, ntcini
+ cova(npvar+(it-1)*ncini+1:it*ncini,npvar+(jt-1)*ncini+1:jt*ncini) = &
+ cova(npvar+(it-1)*ncini+1:it*ncini,npvar+(jt-1)*ncini+1:jt*ncini) + &
+ covar%cortcini(it,jt)*covar%corcini(:,:)
+ end do
end do
+! do i = 1, ntcini*ncini
+! cova(npvar+i,npvar+i) = config%cinierr**2 + cova(npvar+i,npvar+i)
+! end do
endif
! include cross-correlations
diff --git a/source/mod_covar.f90 b/source/mod_covar.f90
index 0120bed..c92725b 100644
--- a/source/mod_covar.f90
+++ b/source/mod_covar.f90
@@ -31,16 +31,24 @@ module mod_covar
implicit none
private
- public :: covar_t, alloc_covar, dealloc_covar
+ public :: covar_t, alloc_covar, dealloc_covar, matrix_oper
type :: covar_t
- real(kind=8), dimension(:), allocatable :: evals ! eigen values of error covariance matrix
- real(kind=8), dimension(:,:), allocatable :: evecs ! eigen vectors of error covariance matrix
- real(kind=8), dimension(:,:), allocatable :: cort ! temporal correlation matrix
- real(kind=8), dimension(:,:), allocatable :: icort ! inverse of temporal correlation matrix
- real(kind=8), dimension(:,:), allocatable :: sqcort ! square root of temporal correlation matrix
- real(kind=8), dimension(:,:), allocatable :: isqcort ! inverse of sqrt of temporal correlation matrix
+ real(kind=8), dimension(:), allocatable :: evals ! eigen values of error covariance matrix
+ real(kind=8), dimension(:,:), allocatable :: evecs ! eigen vectors of error covariance matrix
+ real(kind=8), dimension(:,:), allocatable :: cort ! temporal correlation matrix
+ real(kind=8), dimension(:,:), allocatable :: icort ! inverse of temporal correlation matrix
+ real(kind=8), dimension(:,:), allocatable :: sqcort ! square root of temporal correlation matrix
+ real(kind=8), dimension(:,:), allocatable :: isqcort ! inverse of sqrt of temporal correlation matrix
+ real(kind=8), dimension(:,:), allocatable :: corcini ! spatial correlation matrix for scalars of initial mixing ratio
+ real(kind=8), dimension(:,:), allocatable :: cortcini ! temporal correlation matrix for scalars of initial mixing ratio
+ real(kind=8), dimension(:,:), allocatable :: icortcini ! inverse of temporal correlation matrix of scalars
+ real(kind=8), dimension(:,:), allocatable :: sqcortcini ! square root of temporal correlation matrix of scalars
+ real(kind=8), dimension(:,:), allocatable :: isqcortcini ! inverse of sqrt of temporal correlation matrix of scalars
+ real(kind=8), dimension(:,:), allocatable :: icorcini ! inverse of covariance matrix of scalars
+ real(kind=8), dimension(:,:), allocatable :: sqcorcini ! square root of covariance matrix of scalars
+ real(kind=8), dimension(:,:), allocatable :: isqcorcini ! inverse of sqrt of covariance matrix of scalars
end type
@@ -87,6 +95,46 @@ module mod_covar
write(logid,*) 'ERROR alloc_covar: not enough memory'
stop
endif
+ allocate( covar%corcini(ncini,ncini), stat = ierr )
+ if ( ierr.ne.0 ) then
+ write(logid,*) 'ERROR alloc_covar: not enough memory'
+ stop
+ endif
+ allocate( covar%sqcorcini(ncini,ncini), stat = ierr )
+ if ( ierr.ne.0 ) then
+ write(logid,*) 'ERROR alloc_covar: not enough memory'
+ stop
+ endif
+ allocate( covar%icorcini(ncini,ncini), stat = ierr )
+ if ( ierr.ne.0 ) then
+ write(logid,*) 'ERROR alloc_covar: not enough memory'
+ stop
+ endif
+ allocate( covar%isqcorcini(ncini,ncini), stat = ierr )
+ if ( ierr.ne.0 ) then
+ write(logid,*) 'ERROR alloc_covar: not enough memory'
+ stop
+ endif
+ allocate( covar%cortcini(ntcini,ntcini), stat = ierr )
+ if ( ierr.ne.0 ) then
+ write(logid,*) 'ERROR alloc_covar: not enough memory'
+ stop
+ endif
+ allocate( covar%sqcortcini(ntcini,ntcini), stat = ierr )
+ if ( ierr.ne.0 ) then
+ write(logid,*) 'ERROR alloc_covar: not enough memory'
+ stop
+ endif
+ allocate( covar%icortcini(ntcini,ntcini), stat = ierr )
+ if ( ierr.ne.0 ) then
+ write(logid,*) 'ERROR alloc_covar: not enough memory'
+ stop
+ endif
+ allocate( covar%isqcortcini(ntcini,ntcini), stat = ierr )
+ if ( ierr.ne.0 ) then
+ write(logid,*) 'ERROR alloc_covar: not enough memory'
+ stop
+ endif
end subroutine alloc_covar
@@ -98,16 +146,113 @@ module mod_covar
type (covar_t), intent(in out) :: covar
- if ( allocated(covar%evals) ) deallocate(covar%evals)
- if ( allocated(covar%evecs) ) deallocate(covar%evecs)
- if ( allocated(covar%cort) ) deallocate(covar%cort)
- if ( allocated(covar%icort) ) deallocate(covar%icort)
- if ( allocated(covar%sqcort) ) deallocate(covar%sqcort)
- if ( allocated(covar%isqcort) ) deallocate(covar%isqcort)
+ if ( allocated(covar%evals) ) deallocate(covar%evals)
+ if ( allocated(covar%evecs) ) deallocate(covar%evecs)
+ if ( allocated(covar%cort) ) deallocate(covar%cort)
+ if ( allocated(covar%icort) ) deallocate(covar%icort)
+ if ( allocated(covar%sqcort) ) deallocate(covar%sqcort)
+ if ( allocated(covar%isqcort) ) deallocate(covar%isqcort)
+ if ( allocated(covar%corcini) ) deallocate(covar%corcini)
+ if ( allocated(covar%sqcorcini) ) deallocate(covar%sqcorcini)
+ if ( allocated(covar%icorcini) ) deallocate(covar%icorcini)
+ if ( allocated(covar%isqcorcini) ) deallocate(covar%isqcorcini)
+ if ( allocated(covar%cortcini) ) deallocate(covar%cortcini)
+ if ( allocated(covar%sqcortcini) ) deallocate(covar%sqcortcini)
+ if ( allocated(covar%icortcini) ) deallocate(covar%icortcini)
+ if ( allocated(covar%isqcortcini)) deallocate(covar%isqcortcini)
end subroutine dealloc_covar
! --------------------------------------------------
+ ! matrix_oper
+ ! --------------------------------------------------
+ !
+ !> matrix_oper
+ !!
+ !! Purpose: calculate eigendecomposition of matrix
+ !! and uses this to calculate matrix
+ !! square root and inverse
+ !!
+ ! --------------------------------------------------
+
+ subroutine matrix_oper(ndim, mat, imat, sqmat, isqmat)
+
+ integer, intent(in) :: ndim
+ real(kind=8), dimension(ndim,ndim), intent(in) :: mat
+ real(kind=8), dimension(ndim,ndim), intent(out):: imat, sqmat, isqmat
+
+ 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 :: evals
+ real(kind=8), dimension(:,:), allocatable :: evecs
+ real(kind=8), dimension(ndim) :: tmp
+ integer :: ierr, info, neval, ind
+ integer :: i, j
+
+ ! mat = U*D*U^(T)
+ lwork = 26*ndim
+ liwork = 10*ndim
+
+ allocate( work(lwork), stat=ierr )
+ if ( ierr.ne.0 ) stop 'ERROR mod_covar: not enough memory'
+ allocate( iwork(liwork), stat=ierr )
+ if ( ierr.ne.0 ) stop 'ERROR mod_covar: not enough memory'
+ allocate( isuppz(2*ndim), stat=ierr )
+ if ( ierr.ne.0 ) stop 'ERROR mod_covar: not enough memory'
+ allocate( evecs(ndim,ndim), stat=ierr )
+ if ( ierr.ne.0 ) stop 'ERROR mod_covar: not enough memory'
+ allocate( evals(ndim), stat=ierr )
+ if ( ierr.ne.0 ) stop 'ERROR mod_covar: not enough memory'
+ allocate( zwork(ndim,ndim), stat=ierr )
+ if ( ierr.ne.0 ) stop 'ERROR mod_covar: not enough memory'
+
+ zwork = mat
+ call dsyevr('V','A','L',ndim,zwork,ndim,1.,1.,1,1,0.,neval,evals,evecs,ndim,isuppz,work,lwork,iwork,liwork,info)
+ if ( info.ne.0 ) then
+ write(logid,*) 'ERROR error_cov: eigenvalue decomposition'
+ stop
+ endif
+
+ ! inverse of matrix
+ imat(:,:) = 0d0
+ do i = 1, ndim
+ do j = 1, ndim
+ tmp(j) = evecs(i,j)*1d0/evals(j)
+ end do
+ imat(i,:) = matmul(tmp,transpose(evecs))
+ end do
+
+ ! square root of matrix
+ sqmat(:,:) = 0d0
+ do i = 1, ndim
+ do j = 1, ndim
+ tmp(j) = evecs(i,j)*sqrt(evals(j))
+ end do
+ sqmat(i,:) = matmul(tmp,transpose(evecs))
+ end do
+
+ ! inverse of square root of matrix
+ isqmat(:,:) = 0d0
+ do i = 1, ndim
+ do j = 1, ndim
+ tmp(j) = evecs(i,j)*1d0/sqrt(evals(j))
+ end do
+ isqmat(i,:) = matmul(tmp,transpose(evecs))
+ end do
+
+ deallocate(work)
+ deallocate(zwork)
+ deallocate(iwork)
+ deallocate(isuppz)
+ deallocate(evals)
+ deallocate(evecs)
+
+ end subroutine matrix_oper
+
+ ! --------------------------------------------------
end module mod_covar
diff --git a/source/mod_perturb.f90 b/source/mod_perturb.f90
index 61e8834..301a12b 100644
--- a/source/mod_perturb.f90
+++ b/source/mod_perturb.f90
@@ -66,6 +66,7 @@ module mod_perturb
real(kind=8), dimension(neig) :: tmp
real(kind=8), dimension(ndvar) :: tmp1
+ real(kind=8), dimension(ntcini*ncini) :: tmp2
real(kind=8), dimension(nvar) :: perturb
real, dimension(nvar) :: rvec
integer :: seed
@@ -92,7 +93,19 @@ module mod_perturb
end do
end do
if ( config%opt_cini ) then
- perturb(npvar+1:nvar) = config%cinierr * rvec(npvar+1:nvar)
+ do jt = 1, ntcini
+ do k = 1, ncini
+ do j = 1, ntcini
+ tmp2((j-1)*ncini+1:j*ncini) = covar%sqcortcini(jt,j)*covar%sqcorcini(k,:)
+ end do
+ perturb(npvar+(jt-1)*ncini+k) = dot_product(tmp2,dble(rvec(npvar+1:nvar)))
+ end do
+ end do
+! if ( config%lognormal ) then
+! perturb(npvar+1:nvar) = log(1.+config%cinierr) * rvec(npvar+1:nvar)
+! else
+! perturb(npvar+1:nvar) = config%cinierr * rvec(npvar+1:nvar)
+! endif
endif
! add perturbation to state vector
diff --git a/source/mod_save.f90 b/source/mod_save.f90
index a4e6a5d..71c67c9 100644
--- a/source/mod_save.f90
+++ b/source/mod_save.f90
@@ -199,17 +199,19 @@ module mod_save
integer :: lon_varid, lat_varid, time_varid, ndt_varid
integer :: pri_varid, pos_varid, epri_varid, epos_varid, xpos_varid
integer :: fff_varid, ocn_varid
+ integer :: mcpri_varid, mcpos_varid
integer :: ix, jy, n, nt
integer :: flag, ni, ind
integer :: ix1, ix2, jy1, jy2
real(kind=8), dimension(ntstate) :: timeout
real, dimension(ndt) :: ndtout
- real, dimension(nxregrid,nyregrid,ntstate) :: fpri, fpos, epri, epos
+ real, dimension(nxregrid,nyregrid,ntstate) :: fpri, fpos, epri, epos, mcpri, mcpos
real, dimension(ntstate) :: sumpri, sumpos
- real :: area, landarea, totpri, totpos, totpri_post, totpos_post
+ real :: area, landarea, totpos
real, dimension(nxregrid,nyregrid,ntstate,ndt) :: xpos
real, dimension(nxregrid,nyregrid,ntstate) :: fff_reg, focn_reg
real, dimension(:,:), allocatable :: datain, dataout
+ real, dimension(nxregrid,nyregrid) :: areagrid, work
! regrid fluxes and errors
! ------------------------
@@ -236,15 +238,28 @@ module mod_save
end do
end do
end do
+ if ( config%run_mode.eq.2 ) then
+ ! monte carlo ensemble run for posterior uncertainty
+ ! save posterior and perturbed prior state vectors with no redistribution
+ mcpri(:,:,:) = fpri(:,:,:)
+ mcpos(:,:,:) = fpos(:,:,:)
+ endif
+
+ ! redistribute NEE flux to land
+ ! -----------------------------
! redistribute so that flux increments on ocean boxes are attributed to
! land grid cells within the same aggregated box
! note: future improvement to define aggregated grid with clear separation along coastlines
! then this step no longer needed and should be more accurate!
- totpri = 0.
- totpos = 0.
- totpri_post = 0.
- totpos_post = 0.
+
+ do jy = 1, nyregrid
+ areagrid(:,jy) = grid_area(reg_lat(jy)+0.5*rdy, rdy, rdx)
+ end do
+ work = sum(fpos(:,:,:),dim=3)/real(ntstate)
+ totpos = sum(work(:,:)*areagrid(:,:))*3600.*24.*365./1.e9
+ print*, 'Post increment before redistribution = ',totpos
+
do n = 1, nbox
sumpri(:) = 0.
sumpos(:) = 0.
@@ -252,7 +267,7 @@ module mod_save
do jy = 1, nyregrid
do ix = 1, nxregrid
if ( nbox_xy(ix,jy).eq.n ) then
- area = grid_area(reg_lat(jy)+0.5*rdy, rdy, rdx)
+ area = areagrid(ix,jy)
if ( lsm(ix,jy).ge.0.5 ) then
landarea = landarea + area
endif
@@ -261,10 +276,6 @@ module mod_save
endif
end do
end do
- ! testing
- if(.not.isnan(sum(sumpri))) totpri = totpri + sum(sumpri)/real(ntstate)
- if(.not.isnan(sum(sumpos))) totpos = totpos + sum(sumpos)/real(ntstate)
- ! end testing
if (landarea > 0) then
sumpri(:) = sumpri(:)/landarea
sumpos(:) = sumpos(:)/landarea
@@ -286,19 +297,9 @@ module mod_save
end do
end do
end do
- ! testing
- do jy = 1, nyregrid
- area = grid_area(reg_lat(jy)+0.5*rdy, rdy, rdx)
- do ix = 1, nxregrid
- totpri_post = totpri_post + sum(fpri(ix,jy,:))*area/real(ntstate)
- totpos_post = totpos_post + sum(fpos(ix,jy,:))*area/real(ntstate)
- enddo
- end do
- ! testing
- print*, 'Prior total before redistribution = ',totpri*3600.*24.*365./1.e9
- print*, 'Prior total after redistribution = ',totpri_post*3600.*24.*365./1.e9
- print*, 'Post total before redistribution = ',totpos*3600.*24.*365./1.e9
- print*, 'Post total after redistribution = ',totpos_post*3600.*24.*365./1.e9
+ work = sum(fpos(:,:,:),dim=3)/real(ntstate)
+ totpos = sum(work(:,:)*areagrid(:,:))*3600.*24.*365./1.e9
+ print*, 'Post increment after redistribution = ',totpos
! add prior NEE
! -------------
@@ -316,6 +317,11 @@ module mod_save
endif
fpri(:,jy,:) = fpri(:,jy,:) + dataout(:,:)/numscale
fpos(:,jy,:) = fpos(:,jy,:) + dataout(:,:)/numscale
+ if ( config%run_mode.eq.2 ) then
+ ! for monte carlo ensemble fluxes
+ mcpri(:,jy,:) = mcpri(:,jy,:) + dataout(:,:)/numscale
+ mcpos(:,jy,:) = mcpos(:,jy,:) + dataout(:,:)/numscale
+ endif
end do
print*, 'completed NEE'
deallocate( datain )
@@ -398,6 +404,12 @@ module mod_save
call check( nf90_put_att(ncid, fff_varid, units, flxunit) )
call check( nf90_def_var(ncid, xpos_name, nf90_real, xpos_dimids, xpos_varid) )
call check( nf90_put_att(ncid, xpos_varid, units, flxunit) )
+ if ( config%run_mode.eq.2 ) then
+ call check( nf90_def_var(ncid, 'mcpri', nf90_real, dimids, mcpri_varid) )
+ call check( nf90_put_att(ncid, mcpri_varid, units, flxunit) )
+ call check( nf90_def_var(ncid, 'mcpos', nf90_real, dimids, mcpos_varid) )
+ call check( nf90_put_att(ncid, mcpos_varid, units, flxunit) )
+ endif
call check( nf90_enddef(ncid) )
@@ -413,6 +425,10 @@ module mod_save
call check( nf90_put_var(ncid, ocn_varid, focn_reg) )
call check( nf90_put_var(ncid, fff_varid, fff_reg) )
call check( nf90_put_var(ncid, xpos_varid, xpos) )
+ if ( config%run_mode.eq.2 ) then
+ call check( nf90_put_var(ncid, mcpri_varid, mcpri) )
+ call check( nf90_put_var(ncid, mcpos_varid, mcpos) )
+ endif
call check( nf90_close(ncid) )
@@ -703,13 +719,16 @@ module mod_save
real(kind=8), dimension(nday) :: time, timeout
real, dimension(nxregrid,nyregrid,nday,ndgrid) :: fpri, fpos, xpos
real, dimension(nxregrid,nyregrid,ndgrid) :: flx_nee
+ real, dimension(ndgrid) :: sumpri, sumpos
+ real :: area, landarea, totpos
+ real, dimension(nxregrid,nyregrid) :: areagrid, work
do nd = 1, ndgrid
hours(nd) = (nd-1)*24/ndgrid
end do
- ! interpolate flux increments and add diurnal cycle
- ! -------------------------------------------------
+ ! interpolate flux increments in time
+ ! -----------------------------------
fpri(:,:,:,:) = 0.
fpos(:,:,:,:) = 0.
@@ -722,8 +741,6 @@ module mod_save
d = d + 1
n = int(floor((jd - juldatei)/real(statres,kind=8))) + 1
n = min(n,ntstate)
- ! flux NEE current day
- flx_nee = fluxes%flxnee_nest(:,:,(d-1)*ndgrid+1:d*ndgrid)
do jy = 1, nyregrid
do ix = 1, nxregrid
if ( nbox_xy(ix,jy).le.0 ) cycle
@@ -738,7 +755,7 @@ module mod_save
endif
! note: ndt is number of sub-daily timesteps in state vector
! find index nt to sub-daily timestep in state vector
- do nn = 1, ndt
+ do nn = 1, ndt
if ( (hh.ge.(nn-1)*24/ndt).and.(hh.lt.nn*24/ndt) ) nt = nn
end do
if ( hh.gt.((nt-1)*24/ndt+24/ndt/2).and.nt.lt.ndt ) then
@@ -746,14 +763,8 @@ module mod_save
fpri(ix,jy,d,nd) = states%px0((n-1)*ndt*nbox+(nt-1)*nbox+nbox_xy(ix,jy)) + &
(states%px0((n-1)*ndt*nbox+nt*nbox+nbox_xy(ix,jy)) - &
states%px0((n-1)*ndt*nbox+(nt-1)*nbox+nbox_xy(ix,jy))) * &
- (hh - (nt-1)*24/ndt - 24/ndt/2)*real(ndt)/24. + &
- flx_nee(ix,jy,nd)
+ (hh - (nt-1)*24/ndt - 24/ndt/2)*real(ndt)/24.
fpos(ix,jy,d,nd) = states%px((n-1)*ndt*nbox+(nt-1)*nbox+nbox_xy(ix,jy)) + &
- (states%px((n-1)*ndt*nbox+nt*nbox+nbox_xy(ix,jy)) - &
- states%px((n-1)*ndt*nbox+(nt-1)*nbox+nbox_xy(ix,jy))) * &
- (hh - (nt-1)*24/ndt - 24/ndt/2)*real(ndt)/24. + &
- flx_nee(ix,jy,nd)
- xpos(ix,jy,d,nd) = states%px((n-1)*ndt*nbox+(nt-1)*nbox+nbox_xy(ix,jy)) + &
(states%px((n-1)*ndt*nbox+nt*nbox+nbox_xy(ix,jy)) - &
states%px((n-1)*ndt*nbox+(nt-1)*nbox+nbox_xy(ix,jy))) * &
(hh - (nt-1)*24/ndt - 24/ndt/2)*real(ndt)/24.
@@ -762,24 +773,15 @@ module mod_save
fpri(ix,jy,d,nd) = states%px0((n-1)*ndt*nbox+(nt-2)*nbox+nbox_xy(ix,jy)) + &
(states%px0((n-1)*ndt*nbox+(nt-1)*nbox+nbox_xy(ix,jy)) - &
states%px0((n-1)*ndt*nbox+(nt-2)*nbox+nbox_xy(ix,jy))) * &
- (hh - (nt-2)*24/ndt - 24/ndt/2)*real(ndt)/24. + &
- flx_nee(ix,jy,nd)
+ (hh - (nt-2)*24/ndt - 24/ndt/2)*real(ndt)/24.
fpos(ix,jy,d,nd) = states%px((n-1)*ndt*nbox+(nt-2)*nbox+nbox_xy(ix,jy)) + &
- (states%px((n-1)*ndt*nbox+(nt-1)*nbox+nbox_xy(ix,jy)) - &
- states%px((n-1)*ndt*nbox+(nt-2)*nbox+nbox_xy(ix,jy))) * &
- (hh - (nt-2)*24/ndt - 24/ndt/2)*real(ndt)/24. + &
- flx_nee(ix,jy,nd)
- xpos(ix,jy,d,nd) = states%px((n-1)*ndt*nbox+(nt-2)*nbox+nbox_xy(ix,jy)) + &
(states%px((n-1)*ndt*nbox+(nt-1)*nbox+nbox_xy(ix,jy)) - &
states%px((n-1)*ndt*nbox+(nt-2)*nbox+nbox_xy(ix,jy))) * &
(hh - (nt-2)*24/ndt - 24/ndt/2)*real(ndt)/24.
else
! take current timestep
- fpri(ix,jy,d,nd) = states%px0((n-1)*ndt*nbox+(nt-1)*nbox+nbox_xy(ix,jy)) + &
- flx_nee(ix,jy,nd)
- fpos(ix,jy,d,nd) = states%px((n-1)*ndt*nbox+(nt-1)*nbox+nbox_xy(ix,jy)) + &
- flx_nee(ix,jy,nd)
- xpos(ix,jy,d,nd) = states%px((n-1)*ndt*nbox+(nt-1)*nbox+nbox_xy(ix,jy))
+ fpri(ix,jy,d,nd) = states%px0((n-1)*ndt*nbox+(nt-1)*nbox+nbox_xy(ix,jy))
+ fpos(ix,jy,d,nd) = states%px((n-1)*ndt*nbox+(nt-1)*nbox+nbox_xy(ix,jy))
endif
end do
end do
@@ -787,6 +789,86 @@ module mod_save
time(d) = jd
jd = jd + 1d0
end do
+ ! xpos is simply the increments
+ xpos(:,:,:,:) = fpos(:,:,:,:)
+
+ ! redistribute NEE flux to land
+ ! -----------------------------
+
+ ! redistribute so that flux increments on ocean boxes are attributed to
+ ! land grid cells within the same aggregated box
+ ! note: future improvement to define aggregated grid with clear separation along coastlines
+ ! then this step no longer needed and should be more accurate
+
+ do jy = 1, nyregrid
+ areagrid(:,jy) = grid_area(reg_lat(jy)+0.5*rdy, rdy, rdx)
+ end do
+ work = sum(sum(fpos(:,:,:,:),dim=4),dim=3)/real(nday*ndgrid)
+ totpos = sum(work(:,:)*areagrid(:,:))*3600.*24.*365./numscale/1.e9
+ print*, 'Post increment before redistribution (hourly) = ',totpos
+
+ jd = juldatei
+ d = 0
+ do while (jd.le.juldatef)
+ d = d + 1
+ do n = 1, nbox
+ sumpri(:) = 0.
+ sumpos(:) = 0.
+ landarea = 0.
+ do jy = 1, nyregrid
+ do ix = 1, nxregrid
+ if ( nbox_xy(ix,jy).eq.n ) then
+ area = grid_area(reg_lat(jy)+0.5*rdy, rdy, rdx)
+ if ( lsm(ix,jy).ge.0.5 ) then
+ landarea = landarea + area
+ endif
+ sumpri(:) = sumpri(:) + fpri(ix,jy,d,:)*area
+ sumpos(:) = sumpos(:) + fpos(ix,jy,d,:)*area
+ endif
+ end do
+ end do
+ if (landarea > 0) then
+ sumpri(:) = sumpri(:)/landarea
+ sumpos(:) = sumpos(:)/landarea
+ else
+ sumpri(:) = 0.
+ sumpos(:) = 0.
+ endif
+ do jy = 1, nyregrid
+ do ix = 1, nxregrid
+ if ( nbox_xy(ix,jy).eq.n .and. lsm(ix,jy).ge.0.5 ) then
+ ! land so assign redistributed NEE increment
+ fpri(ix,jy,d,:) = sumpri(:)
+ fpos(ix,jy,d,:) = sumpos(:)
+ else if ( nbox_xy(ix,jy).eq.n ) then
+ ! ocean so NEE increment zero
+ fpri(ix,jy,d,:) = 0.
+ fpos(ix,jy,d,:) = 0.
+ endif
+ end do
+ end do
+ end do ! nbox
+ jd = jd + 1d0
+ end do
+ work = sum(sum(fpos(:,:,:,:),dim=4),dim=3)/real(nday*ndgrid)
+ totpos = sum(work(:,:)*areagrid(:,:))*3600.*24.*365./numscale/1.e9
+ print*, 'Post increment after redistribution (hourly) = ',totpos
+
+ ! add prior NEE
+ ! -------------
+
+ jd = juldatei
+ d = 0
+ do while ( jd.le.juldatef )
+ d = d + 1
+ n = int(floor((jd - juldatei)/real(statres,kind=8))) + 1
+ n = min(n,ntstate)
+ ! flux NEE current day
+ flx_nee = fluxes%flxnee_nest(:,:,(d-1)*ndgrid+1:d*ndgrid)
+ fpri(:,:,d,:) = fpri(:,:,d,:) + flx_nee(:,:,:)
+ fpos(:,:,d,:) = fpos(:,:,d,:) + flx_nee(:,:,:)
+ jd = jd + 1d0
+ end do
fpri(:,:,:,:) = fpri(:,:,:,:)/numscale
fpos(:,:,:,:) = fpos(:,:,:,:)/numscale
xpos(:,:,:,:) = xpos(:,:,:,:)/numscale
diff --git a/source/mod_settings.f90 b/source/mod_settings.f90
index 082047a..535d8ad 100644
--- a/source/mod_settings.f90
+++ b/source/mod_settings.f90
@@ -162,6 +162,7 @@ module mod_settings
real :: sigma_land ! spatial correlation length on land (km)
real :: sigma_ocean ! spatial correlation length on ocean (km)
real :: sigmatime ! temporal correlation length (days)
+ real :: sigmatime_cini ! temporal correlation length for scalars of initial mixing ratios (days)
real :: globerr ! total uncertainty of domain (Tg/y)
end type config_t
@@ -656,6 +657,7 @@ module mod_settings
config%sigma_land = -999.
config%sigma_ocean = -999.
config%sigmatime = -999.
+ config%sigmatime_cini = -999.
config%globerr = -1.
! open file
@@ -863,6 +865,9 @@ module mod_settings
identifier = "sigmatime:"
call read_content (line, identifier, cc, cn, cl, match)
if ( match ) config%sigmatime = real(cn,kind=4)
+ identifier = "sigmatime_cini:"
+ call read_content (line, identifier, cc, cn, cl, match)
+ if ( match ) config%sigmatime_cini = real(cn,kind=4)
identifier = "globerr:"
call read_content (line, identifier, cc, cn, cl, match)
if ( match ) config%globerr = real(cn,kind=4)
@@ -908,41 +913,45 @@ module mod_settings
! initial checks
if ( config%spa_corr ) then
- if ( config%sigma_land.lt.0 ) then
- print*, 'ERROR: uncorrect value of sigma_land'
+ if ( config%sigma_land.lt.0. ) then
+ print*, 'ERROR: incorrect value of sigma_land'
stop
endif
- if ( config%sigma_ocean.lt.0 ) then
- print*, 'ERROR: uncorrect value of sigma_ocean'
+ if ( config%sigma_ocean.lt.0. ) then
+ print*, 'ERROR: incorrect value of sigma_ocean'
stop
endif
endif
- if ( config%sigmatime.lt.0 ) then
- print*, 'ERROR: uncorrect value of sigmatime'
+ if ( config%sigmatime.lt.0. ) then
+ print*, 'ERROR: incorrect value of sigmatime'
+ stop
+ endif
+ if ( config%opt_cini.and.config%sigmatime_cini.lt.0. ) then
+ print*, 'ERROR: incorrect value of sigmatime_cini'
stop
endif
- if ( config%measerr.lt.0 ) then
- print*, 'ERROR: uncorrect value of measerr'
+ if ( config%measerr.lt.0. ) then
+ print*, 'ERROR: incorrect value of measerr'
stop
endif
- if ( config%cinierr.lt.0 ) then
- print*, 'ERROR: uncorrect value of cinierr'
+ if ( config%cinierr.lt.0. ) then
+ print*, 'ERROR: incorrect value of cinierr'
stop
endif
- if ( config%flxerr.lt.0 ) then
- print*, 'ERROR: uncorrect value of flxerr'
+ if ( config%flxerr.lt.0. ) then
+ print*, 'ERROR: incorrect value of flxerr'
stop
endif
- if ( (config%spec.eq.'co2').and.(config%ffferr.lt.0) ) then
- print*, 'ERROR: uncorrect value of ffferr'
+ if ( (config%spec.eq.'co2').and.(config%ffferr.lt.0.) ) then
+ print*, 'ERROR: incorrect value of ffferr'
stop
endif
- if ( config%flxerr_ll.lt.0 ) then
- print*, 'ERROR: uncorrect value of flxerr_ll'
+ if ( config%flxerr_ll.lt.0. ) then
+ print*, 'ERROR: incorrect value of flxerr_ll'
stop
endif
- if ( config%coeff.lt.0 ) then
- print*, 'ERROR: uncorrect value of coeff'
+ if ( config%coeff.lt.0. ) then
+ print*, 'ERROR: incorrect value of coeff'
stop
endif
if ( config%lognormal ) then
diff --git a/source/mod_transform.f90 b/source/mod_transform.f90
index 0fa1038..fe33875 100644
--- a/source/mod_transform.f90
+++ b/source/mod_transform.f90
@@ -43,7 +43,10 @@ module mod_transform
!> chi2phi
!!
!! Purpose: Converts from chi to physical space.
- !! phi = B^(1/2)chi + p0
+ !! For normal PDF:
+ !! phi = B^(1/2)chi + p0
+ !! For lognormal PDF:
+ !! phi = B^(1/2)chi
!!
! --------------------------------------------------
@@ -57,9 +60,12 @@ module mod_transform
real(kind=8), dimension(neig) :: tmp
real(kind=8), dimension(ndvar) :: tmp1
+ real(kind=8), dimension(ntcini*ncini) :: tmp2
integer :: i, j, k, it, jt
chi2phi(:,1) = 0d0
+
+ ! flux variables
do jt = 1, ntstate
do k = 1, ndvar
do j = 1, neig
@@ -74,15 +80,27 @@ module mod_transform
end do
end do
end do
+
+ ! scalars of initial mixing ratio
if ( config%opt_cini ) then
- if ( config%lognormal ) then
- chi2phi(npvar+1:nvar,1) = chi(npvar+1:nvar,1) * log(1.+config%cinierr)
- else
- chi2phi(npvar+1:nvar,1) = chi(npvar+1:nvar,1) * config%cinierr
- endif
+ do jt = 1, ntcini
+ do k = 1, ncini
+ do j = 1, ntcini
+ tmp2((j-1)*ncini+1:j*ncini) = covar%sqcortcini(jt,j)*covar%sqcorcini(k,:)
+ end do
+ chi2phi(npvar+(jt-1)*ncini+k,1) = dot_product(tmp2,chi(npvar+1:nvar,1))
+ end do
+ end do
+! if ( config%lognormal ) then
+! chi2phi(npvar+1:nvar,1) = chi(npvar+1:nvar,1) * log(1.+config%cinierr)
+! else
+! chi2phi(npvar+1:nvar,1) = chi(npvar+1:nvar,1) * config%cinierr
+! endif
endif
- chi2phi(:,1) = chi2phi(:,1) + dble(states%px0)
+ if ( .not.config%lognormal ) then
+ chi2phi(:,1) = chi2phi(:,1) + dble(states%px0)
+ endif
end function chi2phi
@@ -92,7 +110,10 @@ module mod_transform
!> phi2chi
!!
!! Purpose: Converts from physical to chi space.
- !! chi = B^(-1/2)(p - p0)
+ !! For normal PDF:
+ !! chi = B^(-1/2)(p - p0)
+ !! For lognormal PDF:
+ !! chi = B^(-1/2)p
!!
! --------------------------------------------------
@@ -106,9 +127,12 @@ module mod_transform
real(kind=8), dimension(neig) :: tmp
real(kind=8), dimension(ndvar) :: tmp1
+ real(kind=8), dimension(ntcini*ncini) :: tmp2
integer :: i, j, k, it, jt
phi2chi(:,1) = 0d0
+
+ ! flux variables
do jt = 1, ntstate
do k = 1, ndvar
do j = 1, neig
@@ -118,17 +142,36 @@ module mod_transform
tmp1(i) = dot_product(tmp,covar%evecs(i,:))
end do
do it = 1, ntstate
- phi2chi((jt-1)*ndvar+k,1) = phi2chi((jt-1)*ndvar+k,1) + covar%isqcort(jt,it) * &
- dot_product(tmp1,(phi((it-1)*ndvar+1:it*ndvar,1)-dble(states%px0((it-1)*ndvar+1:it*ndvar))))
+ if ( config%lognormal ) then
+ phi2chi((jt-1)*ndvar+k,1) = phi2chi((jt-1)*ndvar+k,1) + covar%isqcort(jt,it) * &
+ dot_product(tmp1,phi((it-1)*ndvar+1:it*ndvar,1))
+ else
+ phi2chi((jt-1)*ndvar+k,1) = phi2chi((jt-1)*ndvar+k,1) + covar%isqcort(jt,it) * &
+ dot_product(tmp1,(phi((it-1)*ndvar+1:it*ndvar,1)-dble(states%px0((it-1)*ndvar+1:it*ndvar))))
+ endif
end do
end do
end do
+
+ ! scalars of initial mixing ratios
if ( config%opt_cini ) then
- if ( config%lognormal ) then
- phi2chi(npvar+1:nvar,1) = (phi(npvar+1:nvar,1) - states%px0(npvar+1:nvar)) * 1./log(1.+config%cinierr)
- else
- phi2chi(npvar+1:nvar,1) = (phi(npvar+1:nvar,1) - states%px0(npvar+1:nvar)) * 1./config%cinierr
- endif
+ do jt = 1, ntcini
+ do k = 1, ncini
+ do j = 1, ntcini
+ tmp2((j-1)*ncini+1:j*ncini) = covar%isqcortcini(jt,j)*covar%isqcorcini(k,:)
+ end do
+ if ( config%lognormal ) then
+ phi2chi(npvar+(jt-1)*ncini+k,1) = dot_product(tmp2,phi(npvar+1:nvar,1))
+ else
+ phi2chi(npvar+(jt-1)*ncini+k,1) = dot_product(tmp2,(phi(npvar+1:nvar,1)-dble(states%px0(npvar+1:nvar))))
+ endif
+ end do
+ end do
+! if ( config%lognormal ) then
+! phi2chi(npvar+1:nvar,1) = phi(npvar+1:nvar,1) * 1./log(1.+config%cinierr)
+! else
+! phi2chi(npvar+1:nvar,1) = (phi(npvar+1:nvar,1) - states%px0(npvar+1:nvar)) * 1./config%cinierr
+! endif
endif
end function phi2chi
@@ -154,9 +197,12 @@ module mod_transform
real(kind=8), dimension(neig) :: tmp
real(kind=8), dimension(ndvar) :: tmp1
+ real(kind=8), dimension(ntcini*ncini) :: tmp2
integer :: i, j, k, it, jt
phi2chi_ad(:,1) = 0d0
+
+ ! flux variables
do jt = 1, ntstate
do k = 1, ndvar
do j = 1, neig
@@ -171,12 +217,22 @@ module mod_transform
end do
end do
end do
+
+ ! scalars of mixing ratios
if ( config%opt_cini ) then
- if ( config%lognormal ) then
- phi2chi_ad(npvar+1:nvar,1) = phi_ad(npvar+1:nvar,1) * log(1.+config%cinierr)
- else
- phi2chi_ad(npvar+1:nvar,1) = phi_ad(npvar+1:nvar,1) * config%cinierr
- endif
+ do jt = 1, ntcini
+ do k = 1, ncini
+ do j = 1, ntcini
+ tmp2((j-1)*ncini+1:j*ncini) = covar%sqcortcini(jt,j)*covar%sqcorcini(k,:)
+ end do
+ phi2chi_ad(npvar+(jt-1)*ncini+k,1) = dot_product(tmp2,phi_ad(npvar+1:nvar,1))
+ end do
+ end do
+! if ( config%lognormal ) then
+! phi2chi_ad(npvar+1:nvar,1) = phi_ad(npvar+1:nvar,1) * log(1.+config%cinierr)
+! else
+! phi2chi_ad(npvar+1:nvar,1) = phi_ad(npvar+1:nvar,1) * config%cinierr
+! endif
endif
end function phi2chi_ad
@@ -199,10 +255,12 @@ module mod_transform
real(kind=8), dimension(neig) :: tmp
real(kind=8), dimension(ndvar) :: tmp1
+ real(kind=8), dimension(ntcini*ncini) :: tmp2
integer :: i, j, k, it, jt
- ! B^(-1)*x for flux variables
prod_icov(:) = 0.
+
+ ! B^(-1)*x for flux variables
do jt = 1, ntstate
do k = 1, ndvar
do j = 1, neig
@@ -217,13 +275,22 @@ module mod_transform
end do
end do
end do
+
! B^(-1)*x for initial mixing ratio scalars
if ( config%opt_cini ) then
- if ( config%lognormal ) then
- prod_icov(npvar+1:nvar) = x(npvar+1:nvar)* 1./(log(1.+config%cinierr))**2
- else
- prod_icov(npvar+1:nvar) = x(npvar+1:nvar)* 1./config%cinierr**2
- endif
+ do jt = 1, ntcini
+ do k = 1, ncini
+ do j = 1, ntcini
+ tmp2((j-1)*ncini+1:j*ncini) = covar%icortcini(jt,j)*covar%icorcini(k,:)
+ end do
+ prod_icov(npvar+(jt-1)*ncini+k) = dot_product(tmp2,x(npvar+1:nvar))
+ end do
+ end do
+! if ( config%lognormal ) then
+! prod_icov(npvar+1:nvar) = x(npvar+1:nvar)* 1./(log(1.+config%cinierr))**2
+! else
+! prod_icov(npvar+1:nvar) = x(npvar+1:nvar)* 1./config%cinierr**2
+! endif
endif
end function prod_icov
diff --git a/source/simulate.f90 b/source/simulate.f90
index 6c8b253..78d72c7 100644
--- a/source/simulate.f90
+++ b/source/simulate.f90
@@ -335,6 +335,7 @@ subroutine simulate(iter, files, config, fluxes, obs, states, grad_o, cost_o)
if ( .not.allocated(istate) ) allocate( istate(ngrid) )
if ( .not.allocated(istateuni) ) allocate( istateuni(ngrid) )
if ( .not.allocated(rstateuni) ) allocate( rstateuni(ngrid) )
+ istate(:) = 0
do n = 1, ngrid
ni = 0
if ( gtime(n).lt.(states%xtime(ntstate)+real(statres,kind=8)) ) then
@@ -428,6 +429,7 @@ subroutine simulate(iter, files, config, fluxes, obs, states, grad_o, cost_o)
if ( config%method.eq.'congrad'.or.config%method.eq.'m1qn3' ) then
! select state variables for current footprint times
allocate( px(ntstep*ndt*nbox), stat=ierr )
+ px(:) = 0.
if( ierr.ne.0 ) then
write(logid,*) 'ERROR simulate: not enough memory'
stop
@@ -469,6 +471,7 @@ subroutine simulate(iter, files, config, fluxes, obs, states, grad_o, cost_o)
! select state variables for current footprint times and
! average transport operator to state variable temporal resolution
allocate( px(ntstep*nbox), stat=ierr )
+ px(:) = 0.
if( ierr.ne.0 ) then
write(logid,*) 'ERROR simulate: not enough memory'
stop
--
GitLab