Commit 839b26d5 authored by ronesy's avatar ronesy
Browse files

Added tool for creating synthetic data for testing and OSEs

parent 7e640f9d
================================================================
FLEXINVERT SYNTHETIC DATA PREPARATION
================================================================
Description:
Prepares synthetic input data for testing FLEXINVERT
Notes:
1) Output from a forward run of FLEXINVERT (run_mode = 0)
are needed before running synthetic data test.
2) For CO2 only NEE fluxes are perturbed.
3) Fluxes are returned on the same domain as the flexpart
grid_time files (if nested is true this is the nested
domain). Requires that the input fluxes cover the same or
greater domain than the flexpart files.
4) If the inversion domain (covered by state vector) is
smaller than the flexpart domain, only the fluxes inside
the inversion domain are perturbed.
Synthetic Data:
1) Prior fluxes: prepared by randomly perturbing prior fluxes
as defined in SETTINGS_[spec]_files using the uncertainty
statistics given in the prior error covariance matrix B
The original prior fluxes (not perturbed) are considered as
the true fluxes.
2) Observations: prepared by forward simulation of FLEXINVERT
using the original (true) prior fluxes. The observations
are randomly perturbed using the uncertainty statistics of
the observation error covariance matrix R.
Usage:
1) Compile with gfortran using: make
2) Edit SETTINGS_files for:
i) path_output: the following files from the forward run
must be in this directory:
- monitor.txt
- nbox_xy.txt
- hloc_box.txt
- evals.txt
- evecs.txt
- cort.txt
The synthetic observations and perturbed prior fluxes
are written to this directory.
ii) path_prior: the prior fluxes are read from this
directory.
Use the same SETTINGS_config as for the forward run.
3) Edit the bash script: job_prep_syndata.sh
4) Run the bash script: ./job_prep_syndata.sh
!---------------------------------------------------------------------------------------
! PREP_SYNDATA: get_error_cov
!---------------------------------------------------------------------------------------
! 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/>.
!
! Copyright 2017, Rona Thompson
!---------------------------------------------------------------------------------------
!
!> get_error_cov
!!
!! Purpose: Reads the eigen values and vectors as well as temporal covariance
!! matrix from an existing FLEXINVERT run.
!!
!---------------------------------------------------------------------------------------
subroutine get_error_cov(files)
use mod_var
use mod_settings
implicit none
type(files_t), intent(in) :: files
character(len=max_path_len) :: filename
character(len=max_path_len) :: rowfmt
logical :: lexist
integer :: ierr, i
real :: tmp
! get number of eigen-values
! --------------------------
filename = trim(files%path_output)//'evals.txt'
inquire(file=trim(filename),exist=lexist)
if (.not.lexist) then
write(logid,*) 'ERROR: cannot find '//trim(filename)
stop
endif
open(100, file=trim(filename), action='read', status='old')
ierr = 0
neig = 0
do while (ierr.eq.0)
read(100,*,iostat=ierr) tmp
neig = neig + 1
end do
close(100)
neig = neig - 1
print*, 'neig = ',neig
allocate ( evals(neig) )
allocate ( evecs(ndvar,neig) )
! read eigen-values
! -----------------
open(100, file=trim(filename), action='read', status='old')
do i = 1, neig
read(100,*) evals(i)
end do
close(100)
! read eigen-vectors
! ------------------
filename = trim(files%path_output)//'evecs.txt'
inquire(file=trim(filename),exist=lexist)
if (.not.lexist) then
write(logid,*) 'ERROR: cannot find '//trim(filename)
stop
endif
write(rowfmt,'(A,I6,A)') '(',neig,'(E11.4,1X))'
open(100, file=trim(filename), action='read', status='old')
do i = 1, ndvar
read(100,fmt=rowfmt) evecs(i,:)
end do
close(100)
! read temporal correlations
! --------------------------
! state vector times
allocate( stat_time(ntstate) )
do i = 1, ntstate
stat_time(i) = juldatei + (i-1)*statres
end do
allocate ( cort(ntstate,ntstate) )
filename = trim(files%path_output)//'cort.txt'
inquire(file=trim(filename),exist=lexist)
if (.not.lexist) then
write(logid,*) 'ERROR: cannot find '//trim(filename)
stop
endif
write(rowfmt,'(A,I6,A)') '(',ntstate,'(E11.4,1X))'
open(100, file=trim(filename), action='read', status='old')
do i = 1, ntstate
read(100,fmt=rowfmt) cort(i,:)
end do
close(100)
end subroutine get_error_cov
!---------------------------------------------------------------------------------------
! PREP_SYNDATA: get_flux
!---------------------------------------------------------------------------------------
! 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/>.
!
! Copyright 2017, Rona Thompson
!---------------------------------------------------------------------------------------
!
!> get_flux
!!
!! Purpose: Reads the prior fluxes and calculates the prior state vector.
!!
!! Interface:
!!
!! Inputs
!! config - configuration settings data structure
!! files - file data structure
!! flx - fluxes
!! state - prior state vector
!!
!! Externals
!! caldate
!! read_ncdf
!!
!---------------------------------------------------------------------------------------
subroutine get_flux(config, files, flx)
use mod_var
use mod_settings
use mod_dates
use mod_tools
use mod_strings
implicit none
type (config_t), intent(in) :: config
type (files_t), intent(in) :: files
real, dimension(nxgrid,nygrid,nt_flx), intent(in out) :: flx
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 :: nread, num
logical :: lexist
real :: area
real :: frac
real, dimension(:,:,:), allocatable :: tmp
allocate( flx_time(nt_flx) )
! read prior fluxes
! -----------------
! variable names
if ( config%spec.eq.'co2' ) then
if ( config%nested ) then
varname = files%varnest_nee
lonname = files%lonnest_nee
latname = files%latnest_nee
timename = files%timenest_nee
else
varname = files%varname_nee
lonname = files%lonname_nee
latname = files%latname_nee
timename = files%timename_nee
endif
else
if ( config%nested ) then
varname = files%varnest_flx
lonname = files%lonnest_flx
latname = files%latnest_flx
timename = files%timenest_flx
else
varname = files%varname_flx
lonname = files%lonname_flx
latname = files%latname_flx
timename = files%timename_flx
endif
endif
n = 0
jd = juldatei
if ( config%spec.eq.'ghg' ) then
! GHG species
do while ( (jd.lt.juldatef).and.(n.lt.nt_flx) )
n = n + 1
call caldate(jd, jjjjmmdd, hhmiss)
write(adate,'(I4.4)') jjjjmmdd/10000
nread = 1
allocate( tmp(nxgrid,nygrid,nread), stat = ierr )
if ( ierr.ne.0 ) stop 'ERROR get_flux: not enough memory'
if ( config%nested ) then
filename = str_replace(files%filenest_flx, 'YYYY', adate)
else
filename = str_replace(files%filename_flx, 'YYYY', adate)
endif
filename = trim(files%path_prior)//trim(filename)
write(logid,*) 'Reading prior fluxes :'//trim(filename)
inquire(file=trim(filename),exist=lexist)
if (.not.lexist) then
write(logid,*) 'ERROR get_flux: cannot find '//trim(filename)
stop
endif
call read_ncdf(filename, &
varname, &
lonname, &
latname, &
timename,&
llx, lly, &
nxgrid, nygrid, &
jd, nread, num, &
tmp)
! convert to kg/m2/s
tmp = tmp/3600.
! apply numerical scaling
tmp = tmp*numscale
flx(:,:,n) = tmp(:,:,1)
deallocate( tmp )
! flux time stamp is start of time step
flx_time(n) = jd
jd = jd + flxtres
end do
else if ( config%spec.eq.'co2' ) then
! CO2 species
do while ( (jd.le.juldatef).and.(n.lt.nt_flx) )
n = n + 1
call caldate(jd, jjjjmmdd, hhmiss)
write(adate,'(I4.4)') jjjjmmdd/10000
! read NEE fluxes for current day (starting from 00:00)
if ( config%nested ) then
filename = str_replace(files%filenest_nee, 'YYYY', adate)
else
filename = str_replace(files%filename_nee, 'YYYY', adate)
endif
print*, 'get_flux: filename = ',filename
filename = trim(files%path_prior)//trim(filename)
inquire(file=trim(filename),exist=lexist)
if (.not.lexist) then
write(logid,*) 'ERROR get_flux: cannot find '//trim(filename)
stop
endif
allocate( tmp(nxgrid,nygrid,nd_nee) )
call read_ncdf(filename, &
varname, &
lonname, &
latname, &
timename,&
llx, lly, &
nxgrid, nygrid, &
jd, nd_nee, num, &
tmp)
! convert to kg/m2/s
tmp = tmp/3600.
! apply numerical scaling
tmp = tmp*numscale
flx(:,:,(n-1)*nd_nee+1:n*nd_nee) = tmp(:,:,:)
deallocate( tmp )
do i = 1, nd_nee
flx_time((n-1)*nd_nee+i) = jd + dble(i-1)/dble(nd_nee)
end do
jd = jd + 1d0
end do
endif
end subroutine get_flux
!---------------------------------------------------------------------------------------
! PREP_SYNDATA: get_obs
!---------------------------------------------------------------------------------------
! 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/>.
!
! Copyright 2017, Rona Thompson
!---------------------------------------------------------------------------------------
!
!> get_obs
!!
!! Purpose: Reads the mixing ratios from the monitor.txt file of an existing
!! FLEXINVERT run.
!!
!! Interface:
!!
!! Inputs
!! config - configuration settings data structure
!! files - file data structure
!! obs - prior modelled mixing ratios (zero vector on input)
!! obserr - error in the observation space (zero vector on input)
!!
!---------------------------------------------------------------------------------------
subroutine get_obs(config, files, obs, obserr)
use mod_var
use mod_settings
implicit none
type(config_t), intent(in) :: config
type(files_t), intent(in) :: files
real, dimension(maxobs), intent(in out) :: obs
real, dimension(maxobs), intent(in out) :: obserr
character(len=max_path_len) :: filename
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
character(len=3) :: recs
obs(:) = 0.
cini(:) = 0.
bkg(:) = 0.
nee(:) = 0.
fff(:) = 0.
ocn(:) = 0.
bbg(:) = 0.
prior(:) = 0.
ghg(:) = 0.
filename = trim(files%path_output)//'monitor.txt'
inquire(file=trim(filename), exist=lexist)
if ( .not.lexist ) then
write(logid,*) 'ERROR: cannot find '//trim(filename)
stop
endif
ierr = 0
open(100,file=trim(filename),status='old',action='read',iostat=ierr)
read(100,fmt='(A)') dump
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)
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))'
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
if ( ierr.ne.0 ) exit
end do
endif
close(100)
nobs = i - 1
print*, 'nobs: ',nobs
if ( config%spec.eq.'co2' ) then
obs = cini + bkg + nee + fff + ocn + bbg
else
if ( config%offsets ) then
obs = cini + bkg + ghg + prior
else
obs = cini + bkg + ocn + prior
endif
endif
end subroutine get_obs
!---------------------------------------------------------------------------------------
! PREP_SYNDATA: initialize
!---------------------------------------------------------------------------------------
! 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/>.
!
! Copyright 2017, Rona Thompson
!---------------------------------------------------------------------------------------
!
!> initialize
!!
!! Purpose: Initializes variables
!!
!---------------------------------------------------------------------------------------
subroutine initialize(files, config)
use mod_var
use mod_settings
use mod_dates
use mod_tools
use mod_flexpart, only : read_header
implicit none
type (files_t), intent(in out) :: files
type (config_t), intent(in) :: config
character(max_path_len) :: filename
character(max_path_len) :: rowfmt
logical :: lexist
integer :: yyyymm
character(6) :: adate
character(4) :: ayear
integer :: numpoint
integer :: ibdate, ibtime
integer, dimension(maxpoint) :: releases
real :: area
integer :: i, ix, jy
real :: bndx, bndy, delx, dely
integer :: numx, numy, xshift
integer :: ierr
! open log file
open(logid,file=trim(files%path_output)//'prep_syndata.log',status='replace',action='write',iostat=ierr)
! read receptor list file
! -----------------------
filename = trim(files%file_recept)
inquire(file=trim(filename),exist=lexist)
if (.not.lexist) then
write(logid,*) 'ERROR: cannot find '//trim(filename)
stop
endif
write(logid,*) 'Reading receptor list: '//trim(filename)
call read_reclist(filename)
! initialize flexpart variables
! -----------------------------
! read header
yyyymm=config%datei/100
write(adate,'(i6)') yyyymm
lexist = .false.
i = 1
do while ( (.not.lexist).and.(i.le.nrec) )
if ( .not.config%nested ) filename = trim(files%path_flexpart)//trim(recname(i))//'/'//adate//'/header'
if ( config%nested ) filename = trim(files%path_flexpart)//trim(recname(i))//'/'//adate//'/header_nest'
inquire(file=trim(filename),exist=lexist)
i = i + 1
end do
if ( .not.lexist ) then
write(logid,*) 'ERROR initialize: cannot find flexpart header'
stop
endif
write(logid,*) 'Reading flexpart header: '//trim(filename)
call read_header(filename, numpoint, ibdate, ibtime, releases, &
numx, numy, bndx, bndy, delx, dely, xshift)
! global domain variables
nxgrid = numx
nygrid = numy
llx = bndx
lly = bndy
dx = delx
dy = dely
nxshift = xshift
if ( .not.allocated(gbl_lon) ) allocate( gbl_lon(nxgrid) )
if ( .not.allocated(gbl_lat) ) allocate( gbl_lat(nygrid) )
do i = 1, nxgrid
gbl_lon(i) = llx + (i-1)*dx
end do
do i = 1, nygrid
gbl_lat(i) = lly + (i-1)*dy
end do
! initialize date variables
! -------------------------
! julian dates for inversion interval
juldatei = juldate(config%datei, 0)
juldatef = juldate(config%datef, 0)
! number days in inversion interval
nday = int(juldatef - juldatei + 1.)
! number of NEE fields per day
if ( config%spec.eq.'co2' ) then
nd_nee = config%num_nee_day
endif