diff --git a/README.txt b/README.txt
index 83cfac9b172d03a329bd9f1af5faac77c50886fa..04d986c1f181593e55356c907f60bbbdc665f985 100644
--- a/README.txt
+++ b/README.txt
@@ -53,4 +53,10 @@ Contents:
grid_to_ncdf:
- extra tool (not needed for FLEXINVERT) to convert
the binary grid_time files to NetCDF
-
+
+ constraint:
+ - extra tool (not necessarily needed for FLEXINVERT) to
+ constrain fluxes to known values
+ - either constrains fluxes > 0 or adds additional constraint
+ for global inversions with limited observational coverage in
+ some regions
diff --git a/constraint/README_constraint.txt b/constraint/README_constraint.txt
new file mode 100644
index 0000000000000000000000000000000000000000..721733b5d1d220a313b211528a3f56c59a52bfca
--- /dev/null
+++ b/constraint/README_constraint.txt
@@ -0,0 +1,30 @@
+================================================================
+
+ FLEXINVERT CONSTRAINT POST-PROCESSOR
+
+================================================================
+
+Description:
+ Constrains FLEXINVERT posterior fluxes to other known values
+ based on the posterior error covariances.
+ Either constrain values of individual grid boxes to known
+ values (e.g. only positive fluxes i.e. ineq:0) OR constrain
+ domain total to a known value (e.g. well known total global
+ emissions from box model).
+
+Note:
+ Setting "const_out" in FLEXINVERT SETTINGS_config file needs
+ to be true to generate output needed for constraint routine.
+
+Usage:
+ 1) run FLEXINVERT with setting const_out: .true.
+ 2) compile with gfortran using: make
+ 3) edit the SETTINGS_constraint file and choose constraint
+ type and inequality value if applicable
+ 4) if total global constraint is activated, create a file
+ containing global fluxes and associated errors for each
+ inversion time step and specify the file path in
+ SETTINGS_constraint. An example for the format of this file
+ can be found in the total_global.txt file
+ 5) edit and run the bash script: job_constraint.sh
+ (or alternatively for slurm: sbatch_constraint.sh)
diff --git a/constraint/SETTINGS_constraint b/constraint/SETTINGS_constraint
new file mode 100644
index 0000000000000000000000000000000000000000..39aaae5c16c92540640ef94928eff84ed95721bd
--- /dev/null
+++ b/constraint/SETTINGS_constraint
@@ -0,0 +1,12 @@
+(1) directory to run inversion (dirinvert)
+/path/to/inversion/directory/
+(2) directory for output (dirout)
+/path/to/output/directory/
+(3) name of log file (logfile)
+constraint.log
+(4) scale flux (scaleflx)
+.true.
+(5) constraint type: inequality [0] or total global [1] (contype)
+1
+(6) if contype=0, inequality constraint as REAL (ineq); if contype=1, path to file with total global constraint (tgconpath)
+/path/to/flexinvertplus/constraint/total_global.txt
diff --git a/constraint/bayesian.f90 b/constraint/bayesian.f90
new file mode 100644
index 0000000000000000000000000000000000000000..48bf113b2ab690ee9a31bc5f41b6a9b4aa6122ca
--- /dev/null
+++ b/constraint/bayesian.f90
@@ -0,0 +1,97 @@
+!---------------------------------------------------------------------------------------
+! FLEXINVERT : bayesian
+!---------------------------------------------------------------------------------------
+! Copyright 2013, Rona Thompson
+!
+! 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/>.
+!---------------------------------------------------------------------------------------
+!> Bayesian inversion with inequality constraints
+!! Optimal value of x:
+!!
+!! x = x0 + AP^T(PAP^T)^-1(c-Px0)
+!!
+!! where:
+!! A = the posterior error covariance matrix
+!! x0 = the posterior state vector without the inequality constraint
+!! P = the constraint selector matrix
+!! c = the constraint vector
+!!
+!! Ref: Thacker, Ocean Modelling, 2007
+!---------------------------------------------------------------------------------------
+
+subroutine bayesian()
+
+ use mod_var
+
+ implicit none
+
+ integer :: i,j,N,lda,info
+ real, dimension(nineq) :: xerr
+ real, dimension(nvar,nineq) :: tselmat
+ real, dimension(nvar,nineq) :: atp,atpipatpc
+ real, dimension(nineq,nineq) :: patp, patpc
+ real, dimension(nineq,nineq) :: tgerrmat
+ real(kind=8), dimension(nineq,nineq) :: M,ipatpc
+
+ allocate(x(nvar))
+ allocate(covmat(nvar,nvar))
+
+ ! calculate (c-Px0)
+ xerr=ineqvec-matmul(selmat,x0)
+
+ ! calculate PAP^T+C
+ tgerrmat(:,:)=0
+ if (contype.eq.1) then
+ do i=1,nineq
+ tgerrmat(i,i)=tgerrvec(i)
+ end do
+ end if
+ tselmat=transpose(selmat)
+ atp=matmul(cova(1:nvar,1:nvar),tselmat)
+ patp=matmul(selmat,atp)
+ patpc=patp+tgerrmat
+
+ ! calculate (PAP^T+C)^-1
+ ! use Cholesky factorization
+ N=nineq ! order of PAP^T
+ M=real(patpc,kind=8) ! symmetric matrix
+ lda=nineq ! leading dimension of PAP^T
+ call dpotrf('L',N,M,lda,info)
+ if(info.gt.0) write(logid,*) 'Error in calculating lower Cholesky triangle'
+ call dpotri('L',N,M,lda,info)
+ if(info.gt.0) write(logid,*) 'Error in inverting (HBH^T + R)'
+ ! returns lower ('L') triangle of inverse of M
+ ipatpc=M
+ do i=1,N
+ do j=1,i
+ ipatpc(j,i)=M(i,j)
+ end do
+ end do
+
+ ! constrainted state variables
+ atpipatpc=matmul(atp,real(ipatpc,kind=4))
+ !atpipatpc=matmul(atp,ipatpc)
+ !write(logid,*) 'Difference from posterior:'
+ !write(logid,*) matmul(atpipatpc,xerr)
+ x=x0+matmul(atpipatpc,xerr)
+
+ !BP Calculate constrainted posterior covariance matrix A*
+ ! A* = A-KPA
+ ! with K = AP^T(PAP^T+C)^â»1
+ !write(logid,*) 'before covmat calc'
+ !covmat=cova(1:nvar,1:nvar)-matmul(matmul(atpipatpc,selmat),cova(1:nvar,1:nvar))
+
+
+end subroutine bayesian
+
diff --git a/constraint/calcselector.f90 b/constraint/calcselector.f90
new file mode 100644
index 0000000000000000000000000000000000000000..210299dabae633ecf86a4a956466f0baf711da0c
--- /dev/null
+++ b/constraint/calcselector.f90
@@ -0,0 +1,75 @@
+!---------------------------------------------------------------------------------------
+! FLEXINVERT : calcselector
+!---------------------------------------------------------------------------------------
+! Copyright 2013, Rona Thompson
+!
+! 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/>.
+!---------------------------------------------------------------------------------------
+!> Calculates the constraint selector matrix
+!---------------------------------------------------------------------------------------
+
+subroutine calcselector()
+
+ use mod_var
+
+ implicit none
+
+ integer :: i,j
+
+ ! find elements of state vector outside the constraint
+
+ if (contype.eq.0) then
+ write(logid,*) 'Constraining grid cells to ',ineq
+ nineq=0
+ do i=1,nvar
+ if(x0(i).lt.ineq) then
+ nineq=nineq+1
+ endif
+ end do
+ write(logid, *) 'Number of values outside constraint:', nineq
+
+ if(nineq.eq.0) then
+ write(logid,*) 'No values violate the constraint'
+ stop
+ endif
+ else
+ write(logid,*) 'Constraining total global'
+ nineq=nemsteps
+ endif
+
+
+ allocate(selmat(nineq,nvar))
+
+ if (contype.eq.0) then ! define selector matrix if constraining grid cells
+ allocate(ineqvec(nineq))
+ ineqvec(:)=ineq*numscale/outheight ! Unit: kg/s*numscale/outheight
+ selmat(:,:)=0
+
+ j=0
+ do i=1,nvar
+ if(x0(i).lt.ineq) then
+ j=j+1
+ selmat(j,i)=1
+ !print*, 'Select:', j, i
+ endif
+ end do
+ else ! use summation matrix if constraining total global
+ ! ineqvec is already defined in readinput.f90
+ selmat(:,:)=0
+ do i=1,nineq
+ selmat(i,(i-1)*nbox+1:i*nbox)=area_box
+ end do
+ endif
+end subroutine calcselector
+
diff --git a/constraint/gridarea.f90 b/constraint/gridarea.f90
new file mode 100644
index 0000000000000000000000000000000000000000..621f7f4441b59cfec09b442bdc344a0d8d66865a
--- /dev/null
+++ b/constraint/gridarea.f90
@@ -0,0 +1,48 @@
+!---------------------------------------------------------------------------------------
+! FLEXINVERT : gridarea
+!---------------------------------------------------------------------------------------
+! Copyright 2013, Rona Thompson
+!
+! 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/>.
+!---------------------------------------------------------------------------------------
+!> Calculates the area of a grid cell
+!---------------------------------------------------------------------------------------
+
+subroutine gridarea(ylat,dy,dx,area)
+
+ implicit none
+
+ real, intent(in) :: ylat,dy,dx
+ real, parameter :: pi=3.14159265, rearth=6.371e6, pih=pi/180.
+ real :: ylata, ylatp, ylatm, cosfact, cosfactp, cosfactm, hzone
+ real, intent(out) :: area
+
+ ylata=ylat+(0.5)*dy
+ ylatp=ylata+0.5*dy
+ ylatm=ylat
+ if (ylatm.lt.0.and.ylatp.gt.0) then
+ hzone = rearth*pih
+ else
+ cosfact=COS(ylata*pih)*rearth
+ cosfactp=COS(ylatp*pih)*rearth
+ cosfactm=COS(ylatm*pih)*rearth
+ if(cosfactp.lt.cosfactm) then
+ hzone=SQRT(rearth**2-cosfactp**2)-SQRT(rearth**2-cosfactm**2)
+ else
+ hzone=SQRT(rearth**2-cosfactm**2)-SQRT(rearth**2-cosfactp**2)
+ endif
+ endif
+ area=2.*pi*rearth*hzone*dx/360.
+
+end subroutine gridarea
diff --git a/constraint/job_constraint.sh b/constraint/job_constraint.sh
new file mode 100644
index 0000000000000000000000000000000000000000..8ae63f25572bd16b7aa24ac2389a80fa3e1197dc
--- /dev/null
+++ b/constraint/job_constraint.sh
@@ -0,0 +1,6 @@
+#!/bin/bash
+#---------------------------------------------------
+
+file='./SETTINGS_constraint'
+
+./constraint ${file}
diff --git a/constraint/main.f90 b/constraint/main.f90
new file mode 100644
index 0000000000000000000000000000000000000000..3bc73e229a4ef72c50fe8a166a8c5a4b53735525
--- /dev/null
+++ b/constraint/main.f90
@@ -0,0 +1,83 @@
+!---------------------------------------------------------------------------------------
+! FLEXINVERT inequality constraint
+!---------------------------------------------------------------------------------------
+! Copyright 2013, Rona Thompson
+!
+! 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/>.
+!---------------------------------------------------------------------------------------
+!> Guide to use:
+!! - applies an inequality constraint to state variables found using
+!! Bayesian inversion with Gaussian PDF and is equivalent to solving the
+!! problem with a truncated Gaussian PDF
+!! - requires inversion output from FLEXINVERT
+!---------------------------------------------------------------------------------------
+
+program main
+
+ use mod_var
+
+ implicit none
+
+ integer :: ierr
+
+ !---------------------------------------------------------------
+ ! initialization
+ !---------------------------------------------------------------
+
+ call readinitial()
+ open(logid,file=trim(dirout)//'/'//trim(logfile),status='replace',action='write',iostat=ierr)
+ if(ierr.ne.0) then
+ print*, 'Error opening logfile'
+ endif
+
+ !---------------------------------------------------------------
+ ! read inputs
+ !---------------------------------------------------------------
+
+ call readinput()
+
+ !---------------------------------------------------------------
+ ! calculate constraint selector matrix
+ !---------------------------------------------------------------
+
+ call calcselector()
+
+ !---------------------------------------------------------------
+ ! bayesian inversion with inequality constraint
+ !---------------------------------------------------------------
+
+ call bayesian()
+
+ !---------------------------------------------------------------
+ ! map fluxes back to flexpart grid
+ !---------------------------------------------------------------
+
+ ! re-grid posterior emissions on fine grid
+ ! if scaleflx = true : uses spatial distribution of prior fluxes
+
+ allocate(emispost(nlon,nlat,ntime))
+
+ write(logid, *) 'Finished calculation, go to regrid'
+ call regridemissions(x,xb,emispost)
+
+ !---------------------------------------------------------------
+ ! save output
+ !---------------------------------------------------------------
+
+ call save()
+
+ close(logid)
+
+end program main
+
diff --git a/constraint/makefile b/constraint/makefile
new file mode 100644
index 0000000000000000000000000000000000000000..2b2c8411768a1363e9c459e9a64ee39793432104
--- /dev/null
+++ b/constraint/makefile
@@ -0,0 +1,44 @@
+F90 = /apps/sw/ubuntu22.04/gcc-11.4.0/gcc-13.2.0-tkesyophy2o6rjlzknndu3b4oyasvuqm/bin/gfortran
+LIBPATH = /apps/sw/ubuntu22.04/gcc-13.2.0/netcdf-fortran-4.6.1-ubcs4l6pjwpgoeyvz32wpzyzykib6bwm/lib/
+INCPATH = /apps/sw/ubuntu22.04/gcc-13.2.0/netcdf-fortran-4.6.1-ubcs4l6pjwpgoeyvz32wpzyzykib6bwm/include/
+
+LNK = -o
+CMPL = -c
+LIBS = -lnetcdf -lnetcdff -llapack
+
+FFLAGS = -O3 -g -m64 -fbounds-check -mcmodel=medium -fconvert=little-endian -frecord-marker=4 -ffree-form \
+ -fbacktrace
+
+LDFLAGS = $(FFLAGS) -L$(LIBPATH) -I$(INCPATH) $(LIBS)
+
+MAIN = constraint
+
+OBJECTS = mod_var.o main.o readinitial.o readinput.o gridarea.o \
+ calcselector.o bayesian.o regridemissions.o write3Dncdf.o save.o
+
+$(MAIN): $(OBJECTS)
+ $(F90) $(LNK) $(MAIN) $(OBJECTS) $(LIBS)
+mod_var.o: mod_var.f90
+ $(F90) $(LDFLAGS) $(CMPL) mod_var.f90
+main.o: main.f90
+ $(F90) $(LDFLAGS) $(CMPL) main.f90
+readinitial.o: readinitial.f90
+ $(F90) $(LDFLAGS) $(CMPL) readinitial.f90
+readinput.o: readinput.f90
+ $(F90) $(LDFLAGS) $(CMPL) readinput.f90
+gridarea.o: gridarea.f90
+ $(F90) $(LDFLAGS) $(CMPL) gridarea.f90
+calcselector.o: calcselector.f90
+ $(F90) $(LDFLAGS) $(CMPL) calcselector.f90
+bayesian.o: bayesian.f90
+ $(F90) $(LDFLAGS) $(CMPL) bayesian.f90
+regridemissions.o: regridemissions.f90
+ $(F90) $(LDFLAGS) $(CMPL) regridemissions.f90
+write3Dncdf.o: write3Dncdf.f90
+ $(F90) $(LDFLAGS) $(CMPL) write3Dncdf.f90
+save.o: save.f90
+ $(F90) $(LDFLAGS) $(CMPL) save.f90
+
+clean:
+ rm *.o *.mod
+
diff --git a/constraint/mod_var.f90 b/constraint/mod_var.f90
new file mode 100644
index 0000000000000000000000000000000000000000..264dd245920547ab240de75e46c61b75099711ed
--- /dev/null
+++ b/constraint/mod_var.f90
@@ -0,0 +1,55 @@
+!---------------------------------------------------------------------------------------
+! FLEXINVERT : mod_var
+!---------------------------------------------------------------------------------------
+! Copyright 2013, Rona Thompson
+!
+! 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/>.
+!---------------------------------------------------------------------------------------
+
+module mod_var
+
+ implicit none
+
+ ! path constants
+ character(len=250) :: dirinvert,dirout
+ character(len=250) :: logfile
+ character(len=250) :: tgconpath
+
+ ! parameters
+ real, parameter :: numscale=1.e12 ! numeric scaling factor
+ real, parameter :: outheight=100. ! height of first layer in flexpart
+ real, parameter :: dx=1.0 ! lon resolution
+ real, parameter :: dy=1.0 ! lat resolution
+ integer, parameter :: logid=500
+ real :: ineq
+ integer :: nemsteps
+ integer :: contype
+ logical :: scaleflx
+
+ ! variables
+ integer :: nxgrid,nygrid,nbox,nvar
+ integer :: nlon,nlat,ntime
+ integer :: nineq
+ integer, dimension(:,:), allocatable :: nbox_xy
+ real, dimension(:), allocatable :: area_box
+ real, dimension(:), allocatable :: lon,lat,time
+ real, dimension(:,:,:), allocatable :: emisprior,emispost_ana,emispost
+ real, dimension(:,:), allocatable :: cova
+ real, dimension(:,:), allocatable :: covmat
+ real, dimension(:), allocatable :: xb,x0,x
+ real, dimension(:,:), allocatable :: selmat
+ real, dimension(:), allocatable :: ineqvec, tgerrvec
+
+end module mod_var
+
diff --git a/constraint/readinitial.f90 b/constraint/readinitial.f90
new file mode 100644
index 0000000000000000000000000000000000000000..b8403b11a00ad8db6daa55a7bb36aee4ae801b57
--- /dev/null
+++ b/constraint/readinitial.f90
@@ -0,0 +1,49 @@
+!---------------------------------------------------------------------------------------
+! FLEXINVERT : readinitial
+!---------------------------------------------------------------------------------------
+! Copyright 2013, Rona Thompson
+!
+! 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/>.
+!---------------------------------------------------------------------------------------
+!> Reads initialization file SETTINGS_constraint
+!---------------------------------------------------------------------------------------
+
+subroutine readinitial
+
+ use mod_var
+
+ implicit none
+
+ open(100,file='SETTINGS_constraint',action='read',status='old')
+ read(100,*)
+ read(100,fmt='(A250)') dirinvert
+ read(100,*)
+ read(100,fmt='(A250)') dirout
+ read(100,*)
+ read(100,*) logfile
+ read(100,*)
+ read(100,*) scaleflx
+ read(100,*)
+ read(100,*) contype
+ read(100,*)
+ if (contype.eq.0) then
+ read(100,*) ineq
+ else if (contype.eq.1) then
+ read(100,fmt='(A250)') tgconpath
+ else
+ write(logid,*) 'Invalid type of constraint: ', contype
+ write(logid,*) 'Must be either 0 or 1'
+ stop
+ endif
+end subroutine readinitial
diff --git a/constraint/readinput.f90 b/constraint/readinput.f90
new file mode 100644
index 0000000000000000000000000000000000000000..ad5e91ea5afe1f55f8a7e42ba826c2f18ff5b5f7
--- /dev/null
+++ b/constraint/readinput.f90
@@ -0,0 +1,277 @@
+!---------------------------------------------------------------------------------------
+! FLEXINVERT : readinput
+!---------------------------------------------------------------------------------------
+! Copyright 2013, Rona Thompson
+!
+! 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/>.
+!---------------------------------------------------------------------------------------
+!> Reads input data from inversion
+!---------------------------------------------------------------------------------------
+
+subroutine readinput()
+
+ use mod_var
+
+ implicit none
+! character(len=100) :: pn
+! pn= '/jetfs/spack/opt/spack/linux-rhel8-skylake/_avx512/gcc-8.3.1'//&
+! &'/netcdf-fortran-4.5.3-kfd2vkj7m6cremij4nunxwqfcrqjovcf/netcdf.inc'
+ include 'netcdf.inc'
+! include pn
+ integer :: i,j,n,ix,jy,nid,ierr,xid,yid,zid,nx,ny,nz,varid
+ real, dimension(:,:), allocatable :: tmp
+ real :: area_grid
+
+ !---------------------------------------------------------------
+ ! prior flux at fine resolution
+ !---------------------------------------------------------------
+
+ ! open netcdf file
+ write(logid,*) 'Reading analysis.nc'
+ ierr=nf_open(trim(dirinvert)//'/'//'analysis.nc',NF_NOWRITE,nid)
+ if(ierr.ne.0) then
+ write(logid,*) 'Error opening analysis.nc'
+ endif
+
+ ! inquire about variables
+ ierr=nf_inq_dimid(nid,'longitude',xid)
+ if(ierr.ne.0) then
+ write(logid,*) 'Error reading dimension: longitude'
+ endif
+ ierr=nf_inq_dimid(nid,'latitude',yid)
+ if(ierr.ne.0) then
+ write(logid,*) 'Error reading dimension: latitude'
+ endif
+ ierr=nf_inq_dimid(nid,'time',zid)
+ if(ierr.ne.0) then
+ write(logid,*) 'Error reading dimension: time'
+ endif
+
+ ! read dimension sizes
+ ierr=nf_inq_dimlen(nid,xid,nlon)
+ if(ierr.ne.0) then
+ write(logid,*) 'Error reading dimension length'
+ endif
+ ierr=nf_inq_dimlen(nid,yid,nlat)
+ if(ierr.ne.0) then
+ write(logid,*) 'Error reading dimension length'
+ endif
+ ierr=nf_inq_dimlen(nid,zid,ntime)
+ if(ierr.ne.0) then
+ write(logid,*) 'Error reading dimension length'
+ endif
+
+ allocate(lon(nlon))
+ allocate(lat(nlat))
+ allocate(time(ntime))
+ allocate(emisprior(nlon,nlat,ntime))
+ allocate(emispost_ana(nlon,nlat,ntime))
+
+ ! read variables
+ ierr=nf_inq_varid(nid,'longitude',varid)
+ ierr=nf_get_var_real(nid,varid,lon)
+ if(ierr.ne.0) then
+ write(logid,*) 'Error reading longitude'
+ endif
+ ierr=nf_inq_varid(nid,'latitude',varid)
+ ierr=nf_get_var_real(nid,varid,lat)
+ if(ierr.ne.0) then
+ write(logid,*) 'Error reading latitude'
+ endif
+ ierr=nf_inq_varid(nid,'time',varid)
+ ierr=nf_get_var_real(nid,varid,time)
+ if(ierr.ne.0) then
+ write(logid,*) 'Error reading time'
+ endif
+ ierr=nf_inq_varid(nid,'fpri',varid)
+ ierr=nf_get_var_real(nid,varid,emisprior)
+ if(ierr.ne.0) then
+ write(logid,*) 'Error reading emisprior'
+ endif
+ ierr=nf_inq_varid(nid,'fpos',varid)
+ ierr=nf_get_var_real(nid,varid,emispost_ana)
+ if(ierr.ne.0) then
+ write(logid,*) 'Error reading emispost_ana'
+ endif
+
+ ierr=nf_close(nid)
+
+ ! multiply by numscale and outheight
+ emisprior=emisprior*numscale/outheight
+ emispost_ana=emispost_ana*numscale/outheight
+
+ !CGZ Adjust nemsteps to number of timesteps, no need to set in SETTINGS_constraint
+ nemsteps=ntime
+
+ !---------------------------------------------------------------
+ ! grid definition
+ !---------------------------------------------------------------
+
+ ! open netcdf file
+ nx=nlon
+ nxgrid=nx
+ ny=nlat
+ nygrid=ny
+ write(logid,*) 'Reading nbox_xy.txt'
+ !print*, 'nlon:', nlon, ' nx:', nx,'nlat:', nlat, ' ny:', ny
+ allocate(nbox_xy(nx,ny))
+ allocate(tmp(nx,ny))
+
+ open(100,file=trim(dirinvert)//'/nbox_xy.txt',action='read',status='old',iostat=ierr)
+ if(ierr.ne.0) then
+ write(logid,*) 'Error reading nbox_xy.txt'
+ endif
+ do i=1,nx
+ read(100,*) tmp(i,:)
+ !write(logid,*) , tmp(i,:)
+ end do
+ nbox_xy=int(tmp)
+ nbox=maxval(nbox_xy)
+
+ !print*, nbox_xy
+ close(100)
+ deallocate(tmp)
+
+ !---------------------------------------------------------------
+ ! grid box area
+ !---------------------------------------------------------------
+
+ allocate(area_box(nbox))
+
+ open(100,file=trim(dirinvert)//'/area_box.txt',action='read',status='old',iostat=ierr)
+ if(ierr.ne.0) then
+ write(logid,*) 'Error reading area_box.txt'
+ endif
+ do i=1,nbox
+ read(100,*) area_box(i)
+ end do
+ !print*, area_box
+ close(100)
+
+
+ !---------------------------------------------------------------
+ ! posterior error covariance
+ !---------------------------------------------------------------
+ write(logid,*) 'Reading cova.txt or izwork.txt'
+ nvar=nbox*nemsteps
+
+ allocate(cova(nvar,nvar))
+ cova(:,:)=0
+ open(100,file=trim(dirinvert)//'/izwork.txt',action='read',status='old',iostat=ierr)
+ if(ierr.ne.0) then
+ write(logid,*) 'Error reading izwork.txt, reading cova.txt'
+ !CGZ: try reading cova.txt
+ open(100,file=trim(dirinvert)//'/cova.txt',action='read',status='old',iostat=ierr)
+ if(ierr.ne.0) then
+ write(logid,*) 'Error reading cova.txt'
+ endif
+ endif
+ do i=1,nvar
+ read(100,*) cova(i,:)
+ ! print*, cova(i,:)
+ end do
+ close(100)
+
+
+ !---------------------------------------------------------------
+ ! prior state vector
+ !---------------------------------------------------------------
+
+ allocate(xb(nvar))
+
+ write(logid,*) 'Calculating prior from analysis.nc'
+ xb(:) = 0.
+ do n = 1, ntime
+ do ix = 1, nlon
+ do jy = 1, nlat
+ call gridarea(lat(1)+(float(jy)-0.5)*dy,dy,dx,area_grid)
+ if ( nbox_xy(ix,jy).gt.0 ) then ! only boxes over land
+ xb((n-1)*nbox+nbox_xy(ix,jy)) = xb((n-1)*nbox+nbox_xy(ix,jy)) + &
+ emisprior(ix,jy,n) * area_grid
+ endif
+ if (ix.eq.1) then
+ endif
+ end do
+ end do
+ end do
+ do n = 1,nemsteps
+ xb((n-1)*nbox+1:(n)*nbox)=xb((n-1)*nbox+1:(n)*nbox)/area_box(:)
+ end do
+
+ !---------------------------------------------------------------
+ ! posterior state vector
+ !---------------------------------------------------------------
+
+ allocate(x0(nvar))
+
+ !write(logid,*) 'Calculating posterior from analysis.nc'
+ !x0(:) = 0.
+ !do n = 1, ntime
+ ! do ix = 1, nlon
+ ! do jy = 1, nlat
+ ! call gridarea(lat(1)+(float(jy)-0.5)*dy,dy,dx,area_grid)
+ ! if ( nbox_xy(ix,jy).gt.0 ) then ! only boxes over land
+ ! x0((n-1)*nbox+nbox_xy(ix,jy)) = xb((n-1)*nbox+nbox_xy(ix,jy)) + &
+ ! emispost_ana(ix,jy,n) * area_grid
+ ! endif
+ ! if (ix.eq.1) then
+ ! endif
+ ! end do
+ ! end do
+ !end do
+ !do n = 1,nemsteps
+ ! x0((n-1)*nbox+1:(n)*nbox)=x0((n-1)*nbox+1:(n)*nbox)/area_box(:)
+ !end do
+
+ write(logid,*) 'Reading posterior.txt'
+ open(100,file=trim(dirinvert)//'posterior.txt',action='read',status='old',iostat=ierr)
+ if(ierr.ne.0) then
+ write(logid,*) 'Error reading posterior.txt'
+ endif
+ write(logid,*) 'nvar:', nvar
+ !write(logid,*) 'nemsteps:', nemsteps
+ write(logid,*) 'nbox:', nbox
+ do i=1,nvar
+ read(100,*) x0(i)
+ end do
+ close(100)
+ ! multiply by numscale and outheight and add prior since flexinvert optimises offstets
+ x0=x0*numscale/outheight+xb
+
+ !---------------------------------------------------------------
+ ! total global constraint with error
+ !---------------------------------------------------------------
+
+ if (contype.eq.1) then
+ allocate(ineqvec(nemsteps))
+ allocate(tgerrvec(nemsteps))
+ write(logid,*) 'Reading ',trim(tgconpath)
+ open(100,file=trim(tgconpath),action='read',status='old',iostat=ierr)
+ if (ierr.ne.0) then
+ write(logid,*) 'Error reading ',trim(tgconpath)
+ endif
+
+ read(100,*)
+ do i=1,nemsteps
+ ! read total global flux and associated error in unit kg/s*numscale/outheight
+ read(100,*) ineqvec(i), tgerrvec(i)
+ ineqvec(i) = ineqvec(i)*numscale/outheight
+ tgerrvec(i) = (tgerrvec(i)*numscale)**2
+ end do
+ close(100)
+ endif
+
+end subroutine readinput
+
+
diff --git a/constraint/regridemissions.f90 b/constraint/regridemissions.f90
new file mode 100644
index 0000000000000000000000000000000000000000..025d1774aaed7e1b4f39b220981ec2365e155bc6
--- /dev/null
+++ b/constraint/regridemissions.f90
@@ -0,0 +1,98 @@
+!---------------------------------------------------------------------------------------
+! FLEXINVERT : regridemissions
+!---------------------------------------------------------------------------------------
+! Copyright 2013, Rona Thompson
+!
+! 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/>.
+!---------------------------------------------------------------------------------------
+!> Regrids emissions (or emission errors) to the flexpart resolution using the prior
+!! mass distribution within each aggregate box
+!---------------------------------------------------------------------------------------
+
+subroutine regridemissions(datapost,dataprior,dataout)
+
+ use mod_var
+
+ implicit none
+
+ integer :: n,nt,ix,jy
+ real :: outlat0,area_grid
+ real, dimension(nbox*nemsteps) :: dataprior,datapost
+ real, dimension(nxgrid,nygrid,nemsteps) :: massratio,dataout
+
+ ! calculate ratio of mass of each grid cell to mass in box aggregate for the prior emissions
+ ! and use this to scale the posterior emissions in each grid cell
+
+ outlat0=lat(1)
+ dataout(:,:,:)=0.
+
+ ! do n=1,nemsteps
+ ! do jy=1,nygrid
+ ! do ix=1,nxgrid
+ ! if(dataprior((n-1)*nbox+nbox_xy(ix,jy)).gt.0.and.scaleflx) then
+ ! call gridarea(outlat0+(float(jy)-0.5)*dy,dy,dx,area_grid)
+ ! massratio(ix,jy,n)=emisprior(ix,jy,n)*area_grid/&
+ ! &(dataprior((n-1)*nbox+nbox_xy(ix,jy))*area_box(nbox_xy(ix,jy)))
+ ! dataout(ix,jy,n)=massratio(ix,jy,n)*datapost((n-1)*nbox+nbox_xy(ix,jy))*&
+ ! &area_box(nbox_xy(ix,jy))/area_grid
+ ! else
+ ! dataout(ix,jy,n)=datapost((n-1)*nbox+nbox_xy(ix,jy))
+
+ ! endif
+ ! end do
+ ! end do
+ ! end do
+
+ !cgz: problems with code above! datapost was only written for 1:nvar, try alternative without scaling function
+ write(logid, *) 'Start regridding over ', nemsteps, nygrid, nxgrid
+ if (.not.scaleflx) then
+ do n=1,nemsteps
+ do jy=1,nygrid
+ do ix=1,nxgrid
+ if(nbox_xy(ix,jy).gt.0.)then
+ ! write(logid,*) 'Use value at:', nbox_xy(ix,jy), 'for ', ix, jy
+ ! dataout(ix,jy,n)=datapost(nbox_xy(ix,jy)) !one time step only
+ dataout(ix,jy,n)=datapost((n-1)*nbox+nbox_xy(ix,jy))
+ endif
+ end do
+ end do
+ end do
+ else
+ write(logid,*) 'Scaling testing!'
+ do n=1,nemsteps
+ do jy=1,nygrid
+ do ix=1,nxgrid
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!if eingefügt!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ if(nbox_xy(ix,jy).gt.0.)then
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ if(dataprior((n-1)*nbox+nbox_xy(ix,jy)).gt.0.)then
+ call gridarea(outlat0+(float(jy)-0.5)*dy,dy,dx,area_grid)
+ massratio(ix,jy,n)=emisprior(ix,jy,n)*area_grid/&
+ &(dataprior((n-1)*nbox+nbox_xy(ix,jy))*area_box(nbox_xy(ix,jy)))
+ dataout(ix,jy,n)=massratio(ix,jy,n)*datapost((n-1)*nbox+nbox_xy(ix,jy))*&
+ &area_box(nbox_xy(ix,jy))/area_grid
+ !write(logid, *) 'Posterior calculation:', area_grid, emisprior(ix,jy,n),dataprior((n-1)*nbox+nbox_xy(ix,jy))&
+ !,area_box(nbox_xy(ix,jy)), &
+ !massratio(ix,jy,n),datapost((n-1)*nbox+nbox_xy(ix,jy)), dataout(ix,jy,n)
+ else
+ dataout(ix,jy,n)=datapost((n-1)*nbox+nbox_xy(ix,jy))
+ endif
+ endif
+ end do
+ end do
+ end do
+ endif
+
+end subroutine regridemissions
+
diff --git a/constraint/save.f90 b/constraint/save.f90
new file mode 100644
index 0000000000000000000000000000000000000000..eec47ec4e16dd53e7f76ca34258625e72d565aaf
--- /dev/null
+++ b/constraint/save.f90
@@ -0,0 +1,64 @@
+!---------------------------------------------------------------------------------------
+! FLEXINVERT : save
+!---------------------------------------------------------------------------------------
+! Copyright 2013, Rona Thompson
+!
+! 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/>.
+!---------------------------------------------------------------------------------------
+!> Save outputs to file
+!---------------------------------------------------------------------------------------
+
+subroutine save()
+
+ use mod_var
+
+ implicit none
+
+ character(len=100) :: fileout,varname,varunit,comment,rowfmt
+ character(len=20) :: dname1,dname2,dname3
+ integer :: n,ierr,i
+
+ ! posterior state vector on coarse grid (x)
+ open(100,file=trim(dirinvert)//'constraint.txt',status='replace',action='write',iostat=ierr)
+ do n=1,nemsteps*nbox
+ write(100,*) x(n)*outheight/numscale
+ end do
+ close(100)
+
+ ! posterior error covariance matrix (C)
+ !open(100,file=trim(dirinvert)//'constraint_cova.txt',status='replace',action='write',iostat=ierr)
+ !write(rowfmt,'(A,I6,A)') '(',nvar,'(E11.4,1X))'
+ !do i = 1, nvar
+ ! write(100,rowfmt) covmat(i,:)
+ !end do
+ !close(100)
+
+ ! gridded posterior fluxes
+
+ emispost=emispost*outheight/numscale
+ !write(logid,*) 'Posterior constrained:', emispost
+
+ fileout='constraint.nc'
+ varname='emis_post'
+ varunit='kg/m2/s'
+ dname1='longitude'
+ dname2='latitude'
+ dname3='time'
+ comment='posterior emissions with inequality constraint'
+ call write3Dncdf(fileout,varname,varunit,emispost,&
+ &nlon,nlat,ntime,lon,lat,time,&
+ &dname1,dname2,dname3,comment)
+
+end subroutine save
+
diff --git a/constraint/sbatch_prep_flux.sh b/constraint/sbatch_prep_flux.sh
new file mode 100644
index 0000000000000000000000000000000000000000..847057eb24b6c8369f5a9e2d1969dfb63dfc1360
--- /dev/null
+++ b/constraint/sbatch_prep_flux.sh
@@ -0,0 +1,13 @@
+#!/bin/bash
+#---------------------------------------------------
+file='./SETTINGS_constraint'
+#---------------------------------------------------
+
+cat <<EOF > run_job.sh
+#!/bin/bash
+./constraint ${file}
+EOF
+
+sbatch --job-name=constraint --mem-per-cpu=10000 run_job.sh
+
+rm -f run_job.sh
diff --git a/constraint/total_global.txt b/constraint/total_global.txt
new file mode 100644
index 0000000000000000000000000000000000000000..ab7adcb57283b3338b5fdb27fca3ec9da245b5a3
--- /dev/null
+++ b/constraint/total_global.txt
@@ -0,0 +1,5 @@
+total domain emissions (kg/s) total domain error (kg/s) [one line per inversion time step]
+1.0000E-2 1.0000E-3
+1.5000E-2 1.0000E-3
+5.1000E-2 2.0000E-3
+3.2500E-2 3.2500E-3
diff --git a/constraint/write3Dncdf.f90 b/constraint/write3Dncdf.f90
new file mode 100644
index 0000000000000000000000000000000000000000..feb06163c98a4cd76971ea4baab6450c9f74173b
--- /dev/null
+++ b/constraint/write3Dncdf.f90
@@ -0,0 +1,96 @@
+!---------------------------------------------------------------------------------------
+! FLEXINVERT : write3Dncdf
+!---------------------------------------------------------------------------------------
+! Copyright 2013, Rona Thompson
+!
+! 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/>.
+!---------------------------------------------------------------------------------------
+!> Writes a 3D array to a netcdf file
+!---------------------------------------------------------------------------------------
+
+subroutine write3Dncdf(fileout,varname,varunit,var,nd1,nd2,nd3,&
+ &dim1,dim2,dim3,dname1,dname2,dname3,comment)
+
+ use mod_var
+
+ implicit none
+
+ !include '/usr/include/netcdf.inc'
+ include 'netcdf.inc'
+
+ character(*) :: fileout,varname,varunit,comment
+ character(len=350)::filename
+ character(len=20) :: dname1,dname2,dname3
+ integer :: nd1,nd2,nd3,did1,did2,did3
+ real :: dim1(nd1),dim2(nd2),dim3(nd3)
+ integer :: nid,varid,ierr,ind
+ integer :: dimids(3)
+ real :: var(nd1,nd2,nd3)
+
+ ! open netcdf file
+ filename=trim(trim(dirout)//trim(fileout))
+ write(logid,*) 'Writing file: ',filename
+ ierr=nf_create(trim(filename),NF_CLOBBER,nid)
+ if(ierr.ne.0) then
+ write(logid,*) 'Error creating file: ',filename, ierr
+ endif
+
+ ! define dimensions
+ ierr=nf_def_dim(nid,dname1,nd1,did1)
+ ierr=nf_def_dim(nid,dname2,nd2,did2)
+ ierr=nf_def_dim(nid,dname3,nd3,did3)
+ ierr=nf_enddef(nid)
+
+ ! define dimension variables
+ ierr=nf_redef(nid)
+ ierr=nf_def_var(nid,dname1,NF_FLOAT,1,(/did1/),varid)
+ ind=len_trim(dname1)
+ ierr=nf_put_att_text(nid,varid,'title',len_trim(dname1),dname1(1:ind))
+ ierr=nf_enddef(nid)
+ ierr=nf_put_var_real(nid,varid,dim1)
+ ierr=nf_redef(nid)
+ ierr=nf_def_var(nid,dname2,NF_FLOAT,1,(/did2/),varid)
+ ind=len_trim(dname2)
+ ierr=nf_put_att_text(nid,varid,'title',len_trim(dname2),dname2(1:ind))
+ ierr=nf_enddef(nid)
+ ierr=nf_put_var_real(nid,varid,dim2)
+ ierr=nf_redef(nid)
+ ierr=nf_def_var(nid,dname3,NF_FLOAT,1,(/did3/),varid)
+ ind=len_trim(dname3)
+ ierr=nf_put_att_text(nid,varid,'title',len_trim(dname3),dname3(1:ind))
+ ierr=nf_enddef(nid)
+ ierr=nf_put_var_real(nid,varid,dim3)
+
+ ! define dependent variable
+ dimids(:)=(/did1,did2,did3/)
+ ierr=nf_redef(nid)
+ ierr=nf_def_var(nid,varname,NF_FLOAT,3,dimids,varid)
+ ind=len_trim(varname)
+ ierr=nf_put_att_text(nid,varid,'title',len_trim(varname),varname(1:ind))
+ ind=len_trim(varunit)
+ ierr=nf_put_att_text(nid,varid,'unit',len_trim(varunit),varunit(1:ind))
+ ierr=nf_enddef(nid)
+
+ ! insert variable
+ ierr=nf_put_var_real(nid,varid,var)
+
+ ! insert variable attribute
+ ierr=nf_redef(nid)
+ ind=len_trim(comment)
+ ierr=nf_put_att_text(nid,varid,'comment',len_trim(comment),comment(1:ind))
+ ierr=nf_enddef(nid)
+
+ ierr=nf_close(nid)
+
+end subroutine write3Dncdf
diff --git a/settings/SETTINGS_aer_config b/settings/SETTINGS_aer_config
index f0112a00d0a97647631658ff282e36799163a7d0..22b28ca62d2f3f4230c0788d78400708b9d1799f 100644
--- a/settings/SETTINGS_aer_config
+++ b/settings/SETTINGS_aer_config
@@ -67,6 +67,11 @@ restart: 0
# only use for debugging small runs
verbose: .true.
+# Posterior error covariance output
+# only use if needed for constraint post-processing
+# only for analytic method
+const_out: .false.
+
# Use satellite observations (true or false)
satellite: .false.
diff --git a/settings/SETTINGS_ch4sat_config b/settings/SETTINGS_ch4sat_config
index ccc7bfb46854fab19d53decb14056e5eb70870cf..5f08c67a70b153067bbeb4568772a56446d63451 100644
--- a/settings/SETTINGS_ch4sat_config
+++ b/settings/SETTINGS_ch4sat_config
@@ -66,6 +66,11 @@ restart: 0
# only use for debugging small runs
verbose: .false.
+# Posterior error covariance output
+# only use if needed for constraint post-processing
+# only for analytic method
+const_out: .false.
+
# Species ("co2" or "ghg")
spec: ghg
diff --git a/settings/SETTINGS_co2_config b/settings/SETTINGS_co2_config
index a80f8fd11788c9daefb366f3a17db5fefa5442bd..5d7e73b3e9d3dd3a13646025acf62203b9b53c50 100644
--- a/settings/SETTINGS_co2_config
+++ b/settings/SETTINGS_co2_config
@@ -67,6 +67,11 @@ restart: 0
# only use for debugging small runs
verbose: .false.
+# Posterior error covariance output
+# only use if needed for constraint post-processing
+# only for analytic method
+const_out: .false.
+
# Use satellite observations (true or false)
satellite: .false.
diff --git a/settings/SETTINGS_ghg_config b/settings/SETTINGS_ghg_config
index 901327c0a80a4a3ad0291e2e0de9c01e60c074db..ed3eb4077a016c8ba0b34609a1d31c55c5e92bba 100644
--- a/settings/SETTINGS_ghg_config
+++ b/settings/SETTINGS_ghg_config
@@ -66,6 +66,11 @@ restart: 0
# only use for debugging small runs
verbose: .true.
+# Posterior error covariance output
+# only use if needed for constraint post-processing
+# only for analytic method
+const_out: .true.
+
# Use satellite observations (true or false)
satellite: .false.
diff --git a/source/mod_analytic.f90 b/source/mod_analytic.f90
index e79d9f7ef926175a065a91f77d8ea24479639185..4f1d9fa2bd6c4b606208065ee15cf6675e0e60fc 100644
--- a/source/mod_analytic.f90
+++ b/source/mod_analytic.f90
@@ -308,7 +308,7 @@ module mod_analytic
! verbose output
! --------------
- if ( config%verbose ) then
+ if ( config%const_out ) then
filename = trim(files%path_output)//'izwork.txt'
open(100,file=trim(filename),status='replace',action='write',iostat=ierr)
write(rowfmt,'(A,I6,A)') '(',nvar,'(E11.4,1X))'
@@ -533,6 +533,16 @@ module mod_analytic
! verbose output
! --------------
+ if (config%const_out) then
+ filename = trim(files%path_output)//'cova.txt'
+ open(100,file=trim(filename),status='replace',action='write',iostat=ierr)
+ write(rowfmt,'(A,I6,A)') '(',nvar,'(E11.4,1X))'
+ do i = 1, nvar
+ write(100,rowfmt) cova(i,:)
+ end do
+ close(100)
+ endif
+
if ( config%verbose ) then
filename = trim(files%path_output)//'imwork.txt'
open(100,file=trim(filename),status='replace',action='write',iostat=ierr)
diff --git a/source/mod_save.f90 b/source/mod_save.f90
index 9c8b2de1366d3ec33805b732a7de096fd2805413..5e3ca753d7c68aca39295aea9694ee00ced38858 100644
--- a/source/mod_save.f90
+++ b/source/mod_save.f90
@@ -139,6 +139,15 @@ module mod_save
call save_scalar_cini(files, config, states)
endif
+ if (config%const_out) then
+ ! write posterior state vector if constraint output is requested
+ open(100,file=trim(files%path_output)//'posterior.txt',status='replace',action='write',iostat=ierr)
+ do n = 1, nvar
+ write(100,*) states%px(n)/numscale
+ end do
+ close(100)
+ endif
+
end subroutine save_state
! --------------------------------------------------
@@ -450,6 +459,7 @@ module mod_save
integer :: pri_varid, pos_varid, epri_varid, epos_varid, xpos_varid
integer :: mcpri_varid, mcpos_varid
integer :: ix, jy, n, ix0, jy0, nb
+ integer :: ierr
real(kind=8), dimension(ntstate) :: timeout
real, dimension(nxregrid,nyregrid,ntstate) :: fpri, fpos, epri, epos, xpos, focn_reg, mcpri, mcpos
real :: area, totpos
diff --git a/source/mod_settings.f90 b/source/mod_settings.f90
index 4b54d3c3d67bfc518688e8286b88fe0a79ef01cc..d6c74057593bcd71e99fc6c1b3ce30e0c34294e1 100644
--- a/source/mod_settings.f90
+++ b/source/mod_settings.f90
@@ -126,6 +126,7 @@ module mod_settings
logical :: opt_cini ! optimize background mixing ratios (true or false)
logical :: spa_corr ! use spatial correlation (true or false)
logical :: verbose ! verbose output mode (true or false)
+ logical :: const_out ! use output for constraint method (true or false)
integer :: prior_bg ! use prior best guess (0=false, 1=true)
integer :: restart ! restart a run that crashed (0=false, 1=true)
real :: molarmass ! molar mass of species (corresponds to flux files, e.g. C=12, CH4=16)
@@ -642,6 +643,7 @@ module mod_settings
config%opt_cini = .false.
config%spa_corr = .false.
config%verbose = .false.
+ config%const_out = .false.
config%lognormal = .false.
config%satellite = .false.
config%ground = .false.
@@ -732,6 +734,11 @@ module mod_settings
call read_content (line, identifier, cc, cn, cl, match)
if ( match ) config%verbose = cl
+ ! read constraint output setting
+ identifier = "const_out:"
+ call read_content (line, identifier, cc, cn, cl, match)
+ if ( match ) config%const_out = cl
+
! read prior best guess setting
identifier = "prior_bg:"
call read_content (line, identifier, cc, cn, cl, match)