added constraint method with inequality or total global constraint

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 

constraint/README_constraint.txt

+================================================================
+
+ 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)

constraint/SETTINGS_constraint

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

constraint/bayesian.f90

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

constraint/calcselector.f90

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

constraint/gridarea.f90

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

constraint/job_constraint.sh

+#!/bin/bash 
+#---------------------------------------------------
+
+file='./SETTINGS_constraint'
+
+./constraint ${file}

constraint/main.f90

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

constraint/makefile

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

constraint/mod_var.f90

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

constraint/readinitial.f90

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

constraint/readinput.f90

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

constraint/regridemissions.f90

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

constraint/save.f90

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

constraint/sbatch_prep_flux.sh

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

constraint/total_global.txt

+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

constraint/write3Dncdf.f90

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

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.
 

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
 

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.