Commit 53f6d375 authored by ronesy's avatar ronesy
Browse files

Added option to optimize the background mixing ratios for non-CO2 species

parent cb11ccb8
......@@ -125,10 +125,10 @@ subroutine congrad(iter, grad, files, config, fluxes, obs, states, covar)
x0(:,1) = dble(states%px0)
! transform prior state vector chi space
z0 = phi2chi(states, covar, x0)
z0 = phi2chi(config, states, covar, x0)
! transform state vector chi space
zbg = phi2chi(states, covar, xbg)
zbg = phi2chi(config, states, covar, xbg)
! initial lanczos vector is initial gradient in chi space
zgrad0(:,1) = grad(:,1)
......@@ -170,7 +170,7 @@ subroutine congrad(iter, grad, files, config, fluxes, obs, states, covar)
close(100)
else
! transform zw to physical space
zx = chi2phi(states, covar, zw)
zx = chi2phi(config, states, covar, zw)
states%px = real(zx(:,1),kind=4)
cost_o = 0.
grad_o(:) = 0.
......@@ -178,7 +178,7 @@ subroutine congrad(iter, grad, files, config, fluxes, obs, states, covar)
! calculate cost
! Jp = (p - p0)^TB^(-1)(p - p0)
delx = states%px - states%px0
tmp = prod_icov(covar, delx)
tmp = prod_icov(config, covar, delx)
cost_p = dot_product( delx, real(tmp,kind=4) )
! J = 0.5*(Jp + Jo)
cost = 0.5*(cost_p + cost_o)
......@@ -188,7 +188,7 @@ subroutine congrad(iter, grad, files, config, fluxes, obs, states, covar)
! transform gradient to chi space
x_ad(:,1) = dble(grad_o)
z_ad = phi2chi_ad(covar, x_ad )
z_ad = phi2chi_ad(config, covar, x_ad )
! new gradient J' = J'p + J'o
grad(:,1) = zw(:,1) + z_ad(:,1)
......@@ -317,7 +317,7 @@ subroutine congrad(iter, grad, files, config, fluxes, obs, states, covar)
grad(:,1) = grad(:,1) - matmul(zcglwk(:,1:iter), zqg0(1:iter,1))
! transform state vector to physical space
xa = chi2phi(states, covar, za)
xa = chi2phi(config, states, covar, za)
states%px = real(xa(:,1),kind=4)
if (iter.eq.1) go to 10
......@@ -402,6 +402,9 @@ subroutine congrad(iter, grad, files, config, fluxes, obs, states, covar)
end do
end do
end do
if ( config%opt_cini ) then
tmp(npvar+1:nvar) = tmp(npvar+1:nvar) + config%cinierr * pevecs(jk,npvar+1:nvar)
endif
end do
! calculate last part
do i = 1, nvar
......
......@@ -240,10 +240,18 @@ subroutine error_cov(config, files, states, covar, corr, xerr)
end do
end do
do i = 1, nvar
! uncertainty for flux variables
do i = 1, npvar
states%pxerr0(i) = real(sqrt(var(i)),kind=4)
end do
! uncertainty for scalars of initial mixing ratio
if ( config%opt_cini ) then
do i = 1, ntcini*ncini
states%pxerr0(npvar+i) = config%cinierr
end do
endif
! save output
! -----------
......
......@@ -56,11 +56,11 @@ subroutine init_cini(files, obs)
logical :: lexist
integer :: lastjjjjmm, prejjjjmm, projjjjmm
integer :: jjjjmm, jjjjmmdd, hhmiss, mm, eomday
integer :: ix, jy, kz, i
integer :: ix, jy, kz, i, n
real :: conc
real(kind=8) :: jdmid, jdpre, jdlast
real(kind=8), dimension(nobs) :: jdates
real, dimension(nobs) :: cini
real, dimension(nobs,ncini) :: cini
integer, dimension(nobs) :: ind
real, dimension(nxgrid,nygrid,nzgrid) :: gridini
real, dimension(nxgrid,nygrid,nzgrid) :: concini, preconcini, proconcini
......@@ -73,7 +73,7 @@ subroutine init_cini(files, obs)
concini(:,:,:) = 0.
preconcini(:,:,:) = 0.
proconcini(:,:,:) = 0.
cini(:) = 0.
cini(:,:) = 0.
! sort obs chronologically for efficiency
jdates = obs%jdate(:)
......@@ -187,6 +187,12 @@ subroutine init_cini(files, obs)
! calculate contribution of initial concentration to receptor
do kz = 1, nzgrid
do jy = 1, nygrid
n = 0
do while(.true.)
n = n + 1
if ( gbl_lat(jy).lt.cini_lat(n) ) exit
end do
! print*, 'init_cini: gbl_lat(jy), n = ',gbl_lat(jy), n
do ix = 1, nxgrid
! interpolate between concentrations
if( (jdates(i) - trajdays).lt.jdmid ) then
......@@ -196,7 +202,7 @@ subroutine init_cini(files, obs)
conc = concini(ix,jy,kz) + (proconcini(ix,jy,kz) - concini(ix,jy,kz)) * &
real((jdates(i) - trajdays) - jdmid, kind=4)/30.
endif
cini(i) = cini(i) + gridini(ix,jy,kz)*conc
cini(i,n) = cini(i,n) + gridini(ix,jy,kz)*conc
end do
end do
end do
......@@ -211,7 +217,7 @@ subroutine init_cini(files, obs)
! print*, 'init_cini: ind = ',ind
! print*, 'init_cini: cini = ',cini
do i = 1, nobs
obs%cini(ind(i)) = cini(i)
obs%cini(ind(i),:) = cini(i,:)
end do
! print*, 'init_cini: obs%cini = ',obs%cini
......
......@@ -546,19 +546,34 @@ subroutine init_co2(files, config, fluxes, obs, states, covar)
! ------------------
! number of state variables
nvar = ndt*ntstate*nbox
npvar = ndt*ntstate*nbox
ndvar = ndt*nbox
if ( config%opt_cini ) then
nvar = npvar + ntcini*ncini
else
nvar = npvar
endif
write(logid,*) 'Total number of state variables: ',nvar
write(logid,*) 'Number of flux state variables: ',npvar
write(logid,*) 'Number of state variables per flux time step: ',ndvar
write(logid,*) 'Total number of state variables: ',nvar
call alloc_states(states)
allocate( xerr(ndvar) )
! assign prior state vector
states%px0(:) = 0.
if ( config%opt_cini ) states%px0(npvar+1:nvar) = 1.
! assign time stamps to flux time steps
states%xtime(:) = statime
! assign time stamps to mixing ratio scalar time steps
do n = 1, ntcini
states%cinitime(n) = juldatei + dble((n-1)*cinires)
end do
print*, 'init_ghg: cinitime = ',states%cinitime
! assign best guess of state vector
if ( (config%prior_bg.eq.1).and.(config%method.eq.'congrad') ) then
! from previous inversion
......
......@@ -246,7 +246,7 @@ subroutine init_ghg(files, config, fluxes, obs, states, covar)
! set minimum error
do n = 1, ntstate
do nb = 1, nbox
err_box(nb,n) = max(config%flxerr_ll/3600.*numscale,err_box(nb,n))
err_box(nb,n) = max(config%flxerr_ll*numscale/3600.,err_box(nb,n))
end do
end do
......@@ -283,10 +283,16 @@ subroutine init_ghg(files, config, fluxes, obs, states, covar)
! ------------------
! number of state variables
nvar = ntstate*nbox
npvar = ntstate*nbox
ndvar = nbox
write(logid,*) 'Total number of state variables: ',nvar
if ( config%opt_cini ) then
nvar = npvar + ntcini*ncini
else
nvar = npvar
endif
write(logid,*) 'Number of flux state variables: ',npvar
write(logid,*) 'Number of state variables per flux time step: ',ndvar
write(logid,*) 'Total number of state variables: ',nvar
! allocate state and error vectors
call alloc_states(states)
......@@ -299,10 +305,20 @@ subroutine init_ghg(files, config, fluxes, obs, states, covar)
end do
else
! optimize flux offsets
states%px0 = 0.
states%px0(1:npvar) = 0.
endif
! assign prior initial mixing ratio scalars
if ( config%opt_cini ) states%px0(npvar+1:nvar) = 1.
! time stamp of flux variables
states%xtime(:) = fluxes%time
! time stamp of initial mixing ratio scalars
do n = 1, ntcini
states%cinitime(n) = juldatei + dble((n-1)*cinires)
end do
print*, 'init_ghg: cinitime = ',states%cinitime
! average error over time steps
xerr = abs(sum(err_box(:,:),dim=2)/real(ntstate) )
......
......@@ -200,7 +200,10 @@ subroutine initialize(files, config)
! number of NEE fields per day
nd_nee = config%num_nee_day
if( config%spec.eq.'co2' ) then
! number of initial mixing ratio time steps to optimize
ntcini = max(1, nday/cinires)
if ( config%spec.eq.'co2' ) then
! statres is time interval over which state variables are averaged
! statres_hr is the sub-daily time interval for state variables
! ndt is number of time steps per 24h in prior NEE
......@@ -303,6 +306,11 @@ subroutine initialize(files, config)
nbox_xy(:,:) = 0
lsm_box(:,:) = 0
! initial other variables
! -----------------------
maxiter = config%maxiter
write(logid,*) 'nday: ',nday
write(logid,*) 'statres: ',statres
write(logid,*) 'flxtres: ',flxtres
......
......@@ -174,7 +174,7 @@ subroutine m1qn3_interface(iter, grad, files, config, fluxes, obs, states, covar
f = cost
else
! transform x to physical space
px = chi2phi(states, covar, x)
px = chi2phi(config, states, covar, x)
states%px = real(px(:,1),kind=4)
! run model
cost_o = 0.
......@@ -196,7 +196,7 @@ subroutine m1qn3_interface(iter, grad, files, config, fluxes, obs, states, covar
! transform gradient to chi space
phi(:,1) = dble(grad_o)
chi_ad = phi2chi_ad(covar, phi )
chi_ad = phi2chi_ad(config, covar, phi )
! new gradient J' = J'p + J'o
g(:) = x(:) + chi_ad(:,1)
......
......@@ -169,6 +169,13 @@ module mod_analytic
end do
end do
end do
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
end do
endif
! verbose only
! ------------
......@@ -191,6 +198,13 @@ module mod_analytic
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
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
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))'
......@@ -373,6 +387,9 @@ module mod_analytic
end do
end do
end do
if ( config%opt_cini ) then
zwork(npvar+1:nvar,:) = config%cinierr**2 * transpose(dble(hmat(:,npvar+1:nvar)))
endif
! verbose only
! ------------
......@@ -396,6 +413,9 @@ module mod_analytic
matmul(covar%cort(jt,it)*covb,transpose(dble(hmat(:,(it-1)*ndvar+1:it*ndvar))))
end do
end do
if ( config%opt_cini ) then
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)
write(rowfmt,'(A,I6,A)') '(',nobs,'(E11.4,1X))'
......@@ -496,6 +516,13 @@ module mod_analytic
end do
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,:)))
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)
end do
endif
! include cross-correlations
do i = 1, nvar
......@@ -545,8 +572,15 @@ module mod_analytic
real, intent (in out) :: cost_o
integer :: i
! posterior mixing ratio
obs%model = matmul(hmat,states%px)
! posterior mixing ratio contribution from fluxes in domain
obs%model = matmul(hmat(:,1:npvar),states%px(1:npvar))
! posterior initial mixing ratios
if ( config%opt_cini ) then
obs%cinipos(:) = matmul(hmat(:,npvar+1:nvar),states%px(npvar+1:nvar))
else
obs%cinipos(:) = sum(obs%cini(:,:),dim=2)
endif
! cost observation space
! Jo = (y - Hx)^TR^(-1)(y - Hx)
......@@ -554,14 +588,14 @@ module mod_analytic
do i = 1, nobs
if ( config%spec.eq.'co2' ) then
obs%delta(i) = obs%model(i) + obs%nee(i) + obs%fff(i) + obs%ocn(i) + obs%bbg(i) - &
obs%conc(i) + obs%bkg(i) + obs%cini(i)
obs%conc(i) + obs%bkg(i) + obs%cinipos(i)
else
if ( .not.config%offsets ) then
! optimize fluxes themselves
obs%delta(i) = obs%model(i) + obs%ocn(i) - obs%conc(i) + obs%bkg(i) + obs%cini(i)
obs%delta(i) = obs%model(i) + obs%ocn(i) - obs%conc(i) + obs%bkg(i) + obs%cinipos(i)
else
! optimize offsets so account for contribution from best guess fluxes (ghg)
obs%delta(i) = obs%model(i) + obs%ghg(i) - obs%conc(i) + obs%bkg(i) + obs%cini(i)
obs%delta(i) = obs%model(i) + obs%ghg(i) - obs%conc(i) + obs%bkg(i) + obs%cinipos(i)
endif
endif
cost_o = cost_o + obs%delta(i)**2/obs%err(i)**2
......
......@@ -50,7 +50,8 @@ module mod_obs
real, dimension(:), allocatable :: model ! state vector contribution to mixing ratio
real, dimension(:), allocatable :: delta ! model-observation difference
real, dimension(:), allocatable :: bkg ! background concentrations
real, dimension(:), allocatable :: cini ! inital concentrations
real, dimension(:,:), allocatable :: cini ! inital concentrations
real, dimension(:), allocatable :: cinipos ! posterior inital concentrations
real, dimension(:), allocatable :: measerr ! measurement error
real, dimension(:), allocatable :: err ! total observation error
......@@ -108,7 +109,12 @@ module mod_obs
write(logid,*) 'ERROR alloc_obs: not enough memory'
stop
endif
allocate( obs%cini(nobs), stat = ierr )
allocate( obs%cini(nobs,ncini), stat = ierr )
if ( ierr.ne.0 ) then
write(logid,*) 'ERROR alloc_obs: not enough memory'
stop
endif
allocate( obs%cinipos(nobs), stat = ierr )
if ( ierr.ne.0 ) then
write(logid,*) 'ERROR alloc_obs: not enough memory'
stop
......@@ -154,7 +160,7 @@ module mod_obs
obs%prior(:) = 0.
obs%delta(:) = 0.
obs%bkg(:) = 0.
obs%cini(:) = 0.
obs%cini(:,:) = 0.
obs%measerr(:) = 0.
obs%err(:) = 0.
obs%ocn(:) = 0.
......@@ -225,7 +231,7 @@ module mod_obs
obs%jdate(i) = jdate
obs%conc(i) = conc
obs%measerr(i) = err
obs%err(i) = sqrt(err**2 + config%cinierr**2)
obs%err(i) = err
end do
close(200)
......
......@@ -91,6 +91,9 @@ module mod_perturb
end do
end do
end do
if ( config%opt_cini ) then
perturb(npvar+1:nvar) = config%cinierr * rvec(npvar+1:nvar)
endif
! add perturbation to state vector
states%px0(:) = states%px0(:) + real(perturb(:),kind=4)
......
......@@ -35,7 +35,7 @@ module mod_save
implicit none
private
public :: save_state, save_co2, save_ghg, save_obs, save_nee_hourly, check
public :: save_state, save_co2, save_ghg, save_obs, save_nee_hourly, save_scalar_cini, check
contains
......@@ -64,6 +64,7 @@ module mod_save
character(len=max_path_len) :: rowfmt
integer :: jjjjmmdd, hhmiss
integer :: i, ierr
real, dimension(nobs) :: cinipri
filename = trim(files%path_output)//'monitor.txt'
open(100,file=trim(filename),status='replace',action='write',iostat=ierr)
......@@ -72,28 +73,33 @@ module mod_save
stop
endif
cinipri(:) = sum(obs%cini(:,:),dim=2)
if ( config%spec.eq.'co2' ) then
write(100,*) 'rec yyyymmdd hhmmss juldate conc cini bkg nee fff ocn bbg prior post diff error'
write(100,*) 'rec yyyymmdd hhmmss juldate conc cini cinipos bkg nee fff ocn bbg prior post diff error'
do i = 1, nobs
call caldate(obs%jdate(i), jjjjmmdd, hhmiss)
write(rowfmt,'(A,I6,A)') '(A3,1X,I8,1X,I6,1X,F14.6,1X,',11,'(F10.4,1X))'
write(100,fmt=rowfmt) obs%recs(i), jjjjmmdd, hhmiss, obs%jdate(i), obs%conc(i), obs%cini(i), obs%bkg(i), &
write(rowfmt,'(A,I6,A)') '(A3,1X,I8,1X,I6,1X,F14.6,1X,',12,'(F10.4,1X))'
write(100,fmt=rowfmt) obs%recs(i), jjjjmmdd, hhmiss, obs%jdate(i), &
obs%conc(i), cinipri(i), obs%cinipos(i), obs%bkg(i), &
obs%nee(i), obs%fff(i), obs%ocn(i), obs%bbg(i), obs%prior(i), obs%model(i), obs%delta(i), obs%err(i)
end do
else
if ( .not.config%offsets ) then
write(100,*) 'rec yyyymmdd hhmmss juldate conc cini bkg ocn prior post diff error'
write(100,*) 'rec yyyymmdd hhmmss juldate conc cini cinipos bkg ocn prior post diff error'
else
write(100,*) 'rec yyyymmdd hhmmss juldate conc cini bkg ghg prior post diff error'
write(100,*) 'rec yyyymmdd hhmmss juldate conc cini cinipos bkg ghg prior post diff error'
endif
do i = 1, nobs
call caldate(obs%jdate(i), jjjjmmdd, hhmiss)
write(rowfmt,'(A,I6,A)') '(A3,1X,I8,1X,I6,1X,F14.6,1X,',8,'(F10.4,1X))'
write(rowfmt,'(A,I6,A)') '(A3,1X,I8,1X,I6,1X,F14.6,1X,',9,'(F10.4,1X))'
if ( .not.config%offsets ) then
write(100,fmt=rowfmt) obs%recs(i), jjjjmmdd, hhmiss, obs%jdate(i), obs%conc(i), obs%cini(i), obs%bkg(i), &
write(100,fmt=rowfmt) obs%recs(i), jjjjmmdd, hhmiss, obs%jdate(i), &
obs%conc(i), cinipri(i), obs%cinipos(i), obs%bkg(i), &
obs%ocn(i), obs%prior(i), obs%model(i), obs%delta(i), obs%err(i)
else
write(100,fmt=rowfmt) obs%recs(i), jjjjmmdd, hhmiss, obs%jdate(i), obs%conc(i), obs%cini(i), obs%bkg(i), &
write(100,fmt=rowfmt) obs%recs(i), jjjjmmdd, hhmiss, obs%jdate(i), &
obs%conc(i), cinipri(i), obs%cinipos(i), obs%bkg(i), &
obs%ghg(i), obs%prior(i), obs%model(i), obs%delta(i), obs%err(i)
endif
end do
......@@ -131,6 +137,9 @@ module mod_save
else
call save_ghg(files, config, fluxes, states)
endif
if ( config%opt_cini ) then
call save_scalar_cini(files, states)
endif
end subroutine save_state
......@@ -675,6 +684,38 @@ module mod_save
end subroutine save_nee_hourly
! --------------------------------------------------
! save_scalar_cini
! --------------------------------------------------
!> save_scalar_cini
!!
!! Purpose: Saves the scalars of the initial
!! mixing ratios and their uncertainties
!!
! --------------------------------------------------
subroutine save_scalar_cini(files, states)
use mod_settings
use mod_states
use mod_var
type (files_t), intent(in) :: files
type (states_t), intent(in) :: states
character(len=max_path_len) :: rowfmt
integer :: n, ierr
! write to file
open(100,file=trim(files%path_output)//'scalars_cini.txt',status='replace',action='write',iostat=ierr)
write(rowfmt,'(A,I6,A)') '(',4,'(E11.4,1X))'
write(100,*) 'prior, pri_err, post, pos_err'
do n = 1, ntcini*ncini
write(100,rowfmt) states%px0(npvar+n), states%pxerr0(npvar+n), states%px(npvar+n), states%pxerr(npvar+n)
end do
close(100)
end subroutine save_scalar_cini
! --------------------------------------------------
! check
! --------------------------------------------------
......
......@@ -119,8 +119,10 @@ module mod_settings
integer :: datef ! end date (yyyymmdd)
character(len=3) :: spec ! species ('co2' or 'ghg')
character(len=max_name_len) :: method ! inversion method ('analytic' or 'congrad')
integer :: maxiter ! max number of iterations for congrad and m1qn3
logical :: offsets ! optimize offsets (true or false)
logical :: inc_ocean ! include ocean boxes in inversion (true or false)
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)
integer :: prior_bg ! use prior best guess (0=false, 1=true)
......@@ -646,6 +648,7 @@ module mod_settings
! default values for logical variables
config%offsets = .false.
config%inc_ocean = .false.
config%opt_cini = .false.
config%spa_corr = .false.
config%verbose = .false.
......@@ -694,12 +697,18 @@ module mod_settings
identifier = "method:"
call read_content (line, identifier, cc, cn, cl, match)
if ( match ) config%method = cc
identifier = "maxiter:"
call read_content (line, identifier, cc, cn, cl, match)
if ( match ) config%maxiter = int(cn)
identifier = "offsets:"
call read_content (line, identifier, cc, cn, cl, match)
if ( match ) config%offsets = cl
identifier = "inc_ocean:"
call read_content (line, identifier, cc, cn, cl, match)
if ( match ) config%inc_ocean = cl
identifier = "opt_cini:"
call read_content (line, identifier, cc, cn, cl, match)
if ( match ) config%opt_cini = cl
identifier = "spa_corr:"
call read_content (line, identifier, cc, cn, cl, match)
if ( match ) config%spa_corr = cl
......@@ -805,10 +814,16 @@ module mod_settings
end do read_loop
! convert to m
! convert to metres
config%sigma_land = config%sigma_land*1000.
config%sigma_ocean = config%sigma_ocean*1000.
! if optimization of background chosen for CO2 switch it off
if ( config%opt_cini.and.config%spec.eq.'co2' ) then
write(logid,*) 'WARNING read_config_settings: background optimization should not be used for CO2 -> switching off'
config%opt_cini = .false.
endif
end subroutine read_config_settings
! --------------------------------------------------
......
......@@ -36,11 +36,12 @@ module mod_states
type :: states_t
real, dimension(:), allocatable :: px ! state vector
real, dimension(:), allocatable :: px0 ! prior state vector
real, dimension(:), allocatable :: pxerr0 ! prior uncertainty