Commit cba4a794 authored by ronesy's avatar ronesy
Browse files

Update for lognormal pdf in control space and minor bug fixes

parent b28c96f4
......@@ -14,7 +14,7 @@ Contents:
- selects (and optionally averages) observations and prepares
flexpart runs with releases for each observation
- flexpart grid_time output files written per release
- note need FLEXPART-v10.2beta
- note need at least FLEXPART-v10.2
2) prep_fluxes:
- regrids fluxes (fluxes need to be at the same spatial
resolution as the flexpart grid_time files)
......@@ -38,37 +38,3 @@ Contents:
- extra tool (not needed for FLEXINVERT-Plus) to convert
the binary grid_time files to NetCDF
Testing environment:
Basic test data are provided for four example cases, for both
a CH4 and a CO2 inversion and both with and without a nested
domain. The SETTINGS and FLUXES configuration files are
provided for the examples and only need to be editted for the
root directory. Note to run the tests the user has to provide
meteo data for FLEXPART (these are not included).
TEST_INPUT:
1) FLEXPART:
- directory where FLEXPART options files etc are written
- used in "prep_flexpart"
2) OBS:
- directory containing raw observation files
- used in "prep_flexpart"
3) FLUXES:
- directory containing raw flux files
- used in "prep_fluxes"
TEST_OUTPUT:
1) FLEXOUT:
- directory where FLEXPART grid_time and grid_initial
files are written (in step "prep_flexpart")
2) OBS:
- directory where processed observations are written
(in step "prep_flexpart")
3) FLUXES
- directory where processed fluxes are written
(in step "prep_fluxes")
4) RESULTS
- directory where inversion output are written
......@@ -19,8 +19,8 @@ file_recept: /myfile.txt
# Run dates: format yyyymmdd
# datei = start date
# datef = end date
datei: 20140701
datef: 20140831
datei: 20120101
datef: 20120131
# Nested output (true/false)
lnested: .true.
......
......@@ -15,7 +15,7 @@
# ntime = number of time steps
# timeref = ref date if timestamp is julian days (yyyymmdd) otherwise 0
# timestamp = sec or month (if julian days and timeref is not 0 this is ignored)
file_in: /mypath/TEST_INPUT/FLUXES/GHG/CH4_TOTAL_2012_05x05.nc
file_in: ../TEST_CASE_CH4/INPUT/FLUX/CH4_TOTAL_2012_05x05.nc
varname_in: emisch4
lonname_in: longitude
latname_in: latitude
......@@ -38,17 +38,17 @@ timestamp:
# ntime_out = number of time steps
# llx_out = left longitude of domain
# lly_out = lower latitude of domain
file_out: /mypath/TEST_OUTPUT/FLUXES/GHG/CH4_TOTAL_2012_10x10.nc
file_out: ../TEST_CASE_CH4/INPUT/FLUX/CH4_TOTAL_2012_20x20.nc
varname_out: emisch4
varunit: kgC/m2/h
lonname_out: longitude
latname_out: latitude
timename_out: time
year: 2012
nx_out: 360
ny_out: 180
dx_out: 1.0
dy_out: 1.0
nx_out: 180
ny_out: 90
dx_out: 2.0
dy_out: 2.0
ntime_out: 12
llx_out: -180
lly_out: -90
......
......@@ -24,7 +24,7 @@
!!
!---------------------------------------------------------------------------------------
subroutine get_error_cov(files)
subroutine get_error_cov(files, config)
use mod_var
use mod_settings
......@@ -32,6 +32,7 @@ subroutine get_error_cov(files)
implicit none
type(files_t), intent(in) :: files
type(config_t), intent(in) :: config
character(len=max_path_len) :: filename
character(len=max_path_len) :: rowfmt
......@@ -98,6 +99,8 @@ subroutine get_error_cov(files)
stat_time(i) = juldatei + (i-1)*statres
end do
print*, 'stat_time = ',stat_time
allocate ( cort(ntstate,ntstate) )
filename = trim(files%path_output)//'cort.txt'
inquire(file=trim(filename),exist=lexist)
......
......@@ -27,7 +27,6 @@
!! config - configuration settings data structure
!! files - file data structure
!! flx - fluxes
!! state - prior state vector
!!
!! Externals
!! caldate
......@@ -35,7 +34,7 @@
!!
!---------------------------------------------------------------------------------------
subroutine get_flux(config, files, flx)
subroutine get_flux(config, files, flx, state)
use mod_var
use mod_settings
......@@ -48,20 +47,23 @@ subroutine get_flux(config, files, flx)
type (config_t), intent(in) :: config
type (files_t), intent(in) :: files
real, dimension(nxgrid,nygrid,nt_flx), intent(in out) :: flx
real, dimension(nvar), intent(in out) :: state
character(len=max_path_len) :: filename
character(len=max_name_len) :: varname, lonname, latname, timename
character(len=4) :: adate
real(kind=8) :: jd, jdi
integer :: jjjjmmdd, hhmiss
integer :: n, i, ierr
integer :: n, i, ierr, ix, jy, ix1, jy1, nt
integer :: nread, num
logical :: lexist
real :: area
real :: frac
real, dimension(:,:,:), allocatable :: tmp
real, parameter :: smallnum = 1.e-15
allocate( flx_time(nt_flx) )
flx_time(:) = 0.
! read prior fluxes
! -----------------
......@@ -182,6 +184,31 @@ subroutine get_flux(config, files, flx)
endif
! get state variables
! if ( config%lognormal ) then
! ! coordinates to global domain
! ix1 = minloc(gbl_lon,dim=1,mask=abs(gbl_lon-rllx).eq.minval(abs(gbl_lon-rllx)))
! jy1 = minloc(gbl_lat,dim=1,mask=abs(gbl_lat-rlly).eq.minval(abs(gbl_lat-rlly)))
! do nt = 1, ntstate
! do jy = 1, nyregrid
! do ix = 1, nxregrid
! if(nbox_xy(ix,jy).gt.0) then
! area = grid_area(reg_lat(jy)+0.5*rdy, rdy, rdx)
! state((nt-1)*nbox+nbox_xy(ix,jy)) = state((nt-1)*nbox+nbox_xy(ix,jy)) + flx(ix1+ix-1,jy1+jy-1,nt)*area
! endif
! end do
! end do
! end do
! do nt = 1, ntstate
! state((nt-1)*nbox+1:nt*nbox) = state((nt-1)*nbox+1:nt*nbox)/area_box
! end do
! where ( state.le.0 ) state = smallnum
! state = log(state)
! else
state(:) = 0.
! endif
print*, 'flx_time = ',flx_time
end subroutine get_flux
......
......@@ -48,8 +48,8 @@ subroutine get_obs(config, files, obs, obserr)
character(len=max_path_len) :: rowfmt, dump
logical :: lexist
integer :: jjjjmmdd, hhmiss, i, ierr
real :: jdate, conc, post, delta
real, dimension(maxobs) :: cini, bkg, nee, fff, ocn, bbg, prior, ghg
real :: jdate, conc, post, delta, cinipos
real, dimension(maxobs) :: cini, bkg, nee, fff, ocn, prior, ghg
character(len=3) :: recs
obs(:) = 0.
......@@ -58,7 +58,6 @@ subroutine get_obs(config, files, obs, obserr)
nee(:) = 0.
fff(:) = 0.
ocn(:) = 0.
bbg(:) = 0.
prior(:) = 0.
ghg(:) = 0.
......@@ -75,20 +74,15 @@ subroutine get_obs(config, files, obs, obserr)
if ( config%spec.eq.'co2' ) then
write(rowfmt,'(A,I6,A)') '(A3,1X,I8,1X,I6,1X,F14.6,1X,',11,'(F10.4,1X))'
do i = 1, maxobs
read(100,fmt=rowfmt,iostat=ierr) recs, jjjjmmdd, hhmiss, jdate, conc, cini(i), bkg(i), &
nee(i), fff(i), ocn(i), bbg(i), prior(i), post, delta, obserr(i)
read(100,fmt=rowfmt,iostat=ierr) recs, jjjjmmdd, hhmiss, jdate, conc, cini(i), cinipos, bkg(i), &
nee(i), fff(i), ocn(i), prior(i), post, delta, obserr(i)
if ( ierr.ne.0 ) exit
end do
else
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))'
do i = 1, maxobs
if ( config%offsets ) then
read(100,fmt=rowfmt,iostat=ierr) recs, jjjjmmdd, hhmiss, jdate, conc, cini(i), bkg(i), &
ghg(i), prior(i), post, delta, obserr(i)
else
read(100,fmt=rowfmt,iostat=ierr) recs, jjjjmmdd, hhmiss, jdate, conc, cini(i), bkg(i), &
ocn(i), prior(i), post, delta, obserr(i)
endif
read(100,fmt=rowfmt,iostat=ierr) recs, jjjjmmdd, hhmiss, jdate, conc, cini(i), cinipos, bkg(i), &
ghg(i), prior(i), post, delta, obserr(i)
if ( ierr.ne.0 ) exit
end do
endif
......@@ -96,16 +90,30 @@ subroutine get_obs(config, files, obs, obserr)
nobs = i - 1
print*, 'nobs: ',nobs
! check obserr
do i = 1, nobs
if( obserr(i).eq.9.99 .or. obserr(i).eq.0. ) obserr(i) = config%measerr
end do
if ( config%spec.eq.'co2' ) then
obs = cini + bkg + nee + fff + ocn + bbg
obs = cini + bkg + nee + fff + ocn
else
if ( config%offsets ) then
obs = cini + bkg + ghg + prior
if ( config%lognormal ) then
! lognormal
obs = cini + bkg + prior
else
obs = cini + bkg + ocn + prior
! normal
obs = cini + bkg + ghg + prior
endif
endif
open(100,file=trim(files%path_output)//'obs.txt',status='replace',action='write')
do i = 1, nobs
write(100,fmt='(F10.4)') obs(i)
end do
close(100)
end subroutine get_obs
......@@ -113,6 +113,7 @@ subroutine initialize(files, config)
! number days in inversion interval
nday = int(juldatef - juldatei + 1.)
write(logid,*) 'initialize: nday = ',nday
! number of NEE fields per day
if ( config%spec.eq.'co2' ) then
......
......@@ -25,6 +25,7 @@ SRCS = mod_var.f90 \
get_flux.f90 \
get_error_cov.f90 \
get_obs.f90 \
calc_error_cov.f90 \
prep_syndata.f90
OBJECTS = $(SRCS:.f90=.o)
......
......@@ -102,10 +102,11 @@ module mod_perturb
end do
do it = 1, ntstate
perturb((jt-1)*ndvar+k) = perturb((jt-1)*ndvar+k) + sqcort(jt,it) * &
dot_product(tmp1,dble(rvec((it-1)*ndvar+1:it*ndvar)))
dot_product(tmp1,dble(rvec((it-1)*ndvar+1:it*ndvar)))
end do
end do
end do
print*, 'min/max perturbation = ',minval(perturb), maxval(perturb)
! add perturbation to state vector
state(:) = state(:) + real(perturb(:),kind=4)
......@@ -113,15 +114,24 @@ module mod_perturb
! testing
open(100,file=trim(files%path_output)//'pert_state.txt',status='replace',action='write')
do i = 1, nvar
write(100,fmt='(F14.4)') perturb(i)
write(100,fmt='(F14.8)') perturb(i)
end do
close(100)
open(100,file=trim(files%path_output)//'state_pert.txt',status='replace',action='write')
do i = 1, nvar
write(100,fmt='(F14.8)') state(i)
end do
open(100,file=trim(files%path_output)//'rvec.txt',status='replace',action='write')
do i = 1, nvar
write(100,fmt='(F10.8)') rvec(i)
end do
close(100)
if ( config%lognormal ) then
! lognormal
state(:) = exp(state)
endif
end subroutine perturb_flux
! --------------------------------------------------
......@@ -169,6 +179,11 @@ module mod_perturb
write(100,fmt='(F10.4)') perturb(i)
end do
close(100)
open(100,file=trim(files%path_output)//'obs_pert.txt',status='replace',action='write')
do i = 1, nobs
write(100,fmt='(F10.4)') obs(i)
end do
close(100)
end subroutine perturb_obs
......@@ -207,12 +222,11 @@ module mod_perturb
end do
deallocate(seedvec)
! note: normalization not needed with random_normal function!
! normalize random numbers
! not needed with random_normal - RT 12/20
! rvec = rvec*num/sum(rvec)
! mean of zero
! not needed with random_normal - RT 12/20
! rvec = rvec - 1.
end subroutine randnum
......
......@@ -81,6 +81,8 @@ module mod_save
real(kind=8) :: jd
integer :: yyyymmdd, hhmiss
flx_out(:,:,:) = 0.
! coordinates to global domain
ix1 = minloc(gbl_lon,dim=1,mask=abs(gbl_lon-rllx).eq.minval(abs(gbl_lon-rllx)))
ix2 = minloc(gbl_lon,dim=1,mask=abs(gbl_lon-(rurx-rdx)).eq.minval(abs(gbl_lon-(rurx-rdx))))
......@@ -147,43 +149,81 @@ module mod_save
if ( config%spec.eq.'ghg' ) then
! GHG species
do n = 1, nt_flx
jd = juldatei
n = 1
do while ( (jd.lt.juldatef).and.(n.lt.nt_flx) )
nt = minloc(stat_time,dim=1,mask=abs(stat_time-flx_time(n)).eq.minval(abs(stat_time-flx_time(n))))
nt = min(nt,ntstate)
print*, 'stat_time(nt), flx_time(n) = ',stat_time(nt), flx_time(n)
do jy = 1, nyregrid
do ix = 1, nxregrid
! if inc_ocean is true state vector includes all fluxes, if false need to add ocean fluxes
if ( nbox_xy(ix,jy).gt.0 ) then
if ( stat_time(nt).lt.flx_time(n).and.(nt.gt.1) ) then
! interpolate between previous and current timesteps
flx_out(ix,jy,n) = state((nt-2)*nbox+nbox_xy(ix,jy)) + &
(state((nt-1)*nbox+nbox_xy(ix,jy))-state((nt-2)*nbox+nbox_xy(ix,jy))) * &
(flx_time(n) - stat_time(nt))/statres + &
flx(ix1+ix-1,jy1+jy-1,n)
if ( config%lognormal ) then
! lognormal
! flx_out(ix,jy,n) = state((nt-2)*nbox+nbox_xy(ix,jy)) + &
! (state((nt-1)*nbox+nbox_xy(ix,jy))-state((nt-2)*nbox+nbox_xy(ix,jy))) * &
! (flx_time(n) - stat_time(nt))/statres
flx_out(ix,jy,n) = flx(ix1+ix-1,jy1+jy-1,n) * ( state((nt-2)*nbox+nbox_xy(ix,jy)) + &
(state((nt-1)*nbox+nbox_xy(ix,jy))-state((nt-2)*nbox+nbox_xy(ix,jy))) * &
(flx_time(n) - stat_time(nt))/statres )
else
! normal
flx_out(ix,jy,n) = state((nt-2)*nbox+nbox_xy(ix,jy)) + &
(state((nt-1)*nbox+nbox_xy(ix,jy))-state((nt-2)*nbox+nbox_xy(ix,jy))) * &
(flx_time(n) - stat_time(nt))/statres + &
flx(ix1+ix-1,jy1+jy-1,n)
endif
else if ( stat_time(nt).gt.flx_time(n).and.(nt.lt.ntstate) ) then
! interpolate between current and next timesteps
flx_out(ix,jy,n) = state((nt-1)*nbox+nbox_xy(ix,jy)) + &
(state(nt*nbox+nbox_xy(ix,jy))-state((nt-1)*nbox+nbox_xy(ix,jy))) * &
(stat_time(nt) - flx_time(n))/statres + &
flx(ix1+ix-1,jy1+jy-1,n)
if ( config%lognormal ) then
! lognormal
! flx_out(ix,jy,n) = state((nt-1)*nbox+nbox_xy(ix,jy)) + &
! (state(nt*nbox+nbox_xy(ix,jy))-state((nt-1)*nbox+nbox_xy(ix,jy))) * &
! (stat_time(nt) - flx_time(n))/statres
flx_out(ix,jy,n) = flx(ix1+ix-1,jy1+jy-1,n) * ( state((nt-1)*nbox+nbox_xy(ix,jy)) + &
(state(nt*nbox+nbox_xy(ix,jy))-state((nt-1)*nbox+nbox_xy(ix,jy))) * &
(stat_time(nt) - flx_time(n))/statres )
else
! normal
flx_out(ix,jy,n) = state((nt-1)*nbox+nbox_xy(ix,jy)) + &
(state(nt*nbox+nbox_xy(ix,jy))-state((nt-1)*nbox+nbox_xy(ix,jy))) * &
(stat_time(nt) - flx_time(n))/statres + &
flx(ix1+ix-1,jy1+jy-1,n)
endif
else
! take current timestep
flx_out(ix,jy,n) = state((nt-1)*nbox+nbox_xy(ix,jy)) + flx(ix1+ix-1,jy1+jy-1,n)
if ( config%lognormal ) then
! flx_out(ix,jy,n) = state((nt-1)*nbox+nbox_xy(ix,jy))
flx_out(ix,jy,n) = flx(ix1+ix-1,jy1+jy-1,n) * state((nt-1)*nbox+nbox_xy(ix,jy))
else
flx_out(ix,jy,n) = state((nt-1)*nbox+nbox_xy(ix,jy)) + flx(ix1+ix-1,jy1+jy-1,n)
endif
endif
else
! add ocean fluxes if not optimized
flx_out(ix,jy,n) = flx(ix1+ix-1,jy1+jy-1,n)
endif
end do
end do
jd = jd + statres
n = n + 1
end do
endif ! GHG
print*, 'mod_save: n = ',n
! replace fluxes in inversion domain with perturbed fluxes
flx(ix1:ix2,jy1:jy2,:) = flx_out(:,:,:)
flx(ix1:ix2,jy1:jy2,1:n) = flx_out(:,:,1:n)
! undo numerical scaling
flx = flx/numscale
! convert to kg/m2/h
flx = flx*3600.
print*, 'mod_save: min/max flx = ',minval(flx), maxval(flx)
! adjust time to days since 1-Jan-1900
flx_time = flx_time - juldate(19000101, 0) + 1d0
......
......@@ -95,6 +95,8 @@ module mod_settings
integer :: num_nee_day ! number of NEE time steps per day
real :: statres ! number of flux time steps
real :: statres_hr ! sub-daily temporal resolution of state vector (e.g. 6, 12 or 24h)
real :: measerr ! measurement error
logical :: lognormal ! use lognormal distribution in state space (true or false)
end type config_t
......@@ -457,6 +459,9 @@ module mod_settings
logical :: match, cl
character(len=20), dimension(3) :: temp
! logical defaults
config%lognormal = .false.
! open file
open (100, file = trim (filename), status = 'old', iostat=ierr)
if(ierr.gt.0) then
......@@ -504,6 +509,9 @@ module mod_settings
identifier = "offsets:"
call read_content (line, identifier, cc, cn, cl, match)
if ( match ) config%offsets = cl
identifier = "lognormal:"
call read_content (line, identifier, cc, cn, cl, match)
if ( match ) config%lognormal = cl
! read inversion domain settings
identifier = "w_edge_lon:"
......@@ -536,11 +544,18 @@ module mod_settings
identifier = "statres_hr:"
call read_content (line, identifier, cc, cn, cl, match)
if ( match ) config%statres_hr = real(cn,kind=4)
identifier = "measerr:"
call read_content (line, identifier, cc, cn, cl, match)
if ( match ) config%measerr = real(cn)
endif
end do read_loop
if ( config%spec.eq.'co2' ) then
print*, 'WARNING read_config_settings: for CO2 cannot use lognormal -> changing to normal'
config%lognormal = .false.
endif
end subroutine read_config_settings
......
......@@ -31,11 +31,12 @@ module mod_var
integer, parameter :: max_path_len=200 ! max character length of paths and files
integer, parameter :: max_name_len=50 ! max character length of variable names
integer, parameter :: recname_len=3 ! length of receptor names
real, parameter :: numscale=1.e12 ! numeric scaling factor
real, parameter :: numscale=1.e20 ! numeric scaling factor
integer, parameter :: maxpoint=1000 ! max number of releases per receptor per month
integer, parameter :: maxspec=10 ! max number of species in a flexpart run
integer, parameter :: maxlev=50 ! max number of vertical levels in flexpart output
integer, parameter :: maxobs=10000 ! max number of observations
real, parameter :: trunc=1.e-9 ! truncation threshold for eigenvalues of error covariance matrix
! flexpart variables
! ------------------
......
......@@ -76,11 +76,14 @@ program prep_syndata
call initialize(files, config)
! get prior error covariance
call get_error_cov(files)
call get_error_cov(files, config)
! get prior fluxes
allocate( flx(nxgrid,nygrid,nt_flx))
call get_flux(config, files, flx)
allocate( state(nvar) )
flx(:,:,:) = 0.
state(:) = 0.
call get_flux(config, files, flx, state)
! get observations
allocate( obs(maxobs) )
......@@ -91,8 +94,6 @@ program prep_syndata
! ------------
! perturb fluxes
allocate( state(nvar) )
state(:) = 0.
call perturb_flux(config, files, state)
! perturb obs
......
......@@ -26,14 +26,13 @@ datef: 20120131
# Inversion method ('analytic', 'congrad' or 'm1qn3')
method: m1qn3
# Average/select flexpart releases (true or false)
average_fp: .false.
# Number of iterations
# only used if method is 'congrad' or 'm1qn3'
maxiter: 2
# Optimize flux offsets (true or false)
# optimizes the offsets rather than the fluxes themselves (ghg species only)
offsets: .false.
# Include ocean boxes in inversion (true or false)
# currently only for ghg species
inc_ocean: .false.
......@@ -112,8 +111,8 @@ num_nee_day: 8
# used if error in obs input <= zero
measerr: 0.5
# Initial concentration error: unit same as obs
cinierr: 1.0
# Scalar of initial conc error: fraction
cinierr: 0.01
# Prior flux error: fraction
flxerr: 0.25
......
......@@ -26,14 +26,13 @@ datef: 20120131
# Inversion method ('analytic', 'congrad' or 'm1qn3')
method: analytic
# Average/select flexpart releases (true or false)
average_fp: .false.
# Number of iterations
# only used if method is 'congrad' or 'm1qn3'
maxiter: 2
# Optimize flux offsets (true or false)
# optimizes the offsets rather than the fluxes themselves (ghg species only)
offsets: .false.
# Include ocean boxes in inversion (true or false)
# currently only for ghg species
inc_ocean: .false.
......@@ -112,8 +111,8 @@ num_nee_day: 8
# used if error in obs input <= zero