diff --git a/README.txt b/README.txt
index 0132ef372de00c99b0a72185b8de8619086d4583..fe9fe0a6852fb3d015ba5d96de8f2863463642db 100644
--- a/README.txt
+++ b/README.txt
@@ -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
-
-
-
diff --git a/grid_to_ncdf/SETTINGS b/grid_to_ncdf/SETTINGS
index 0020b7786a6d42fa85c5a893138ed2c5f337d460..c4f81b08ad1776246500a72eb8455271f2e9b606 100644
--- a/grid_to_ncdf/SETTINGS
+++ b/grid_to_ncdf/SETTINGS
@@ -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.
diff --git a/prep_fluxes/FLUXES_ghg b/prep_fluxes/FLUXES_ghg
index f85dbb92658ddb787124d453c3a402988d25a60c..879ec1205f8bcfac43dd4a6086a195e7bab0ae24 100644
--- a/prep_fluxes/FLUXES_ghg
+++ b/prep_fluxes/FLUXES_ghg
@@ -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
diff --git a/prep_regions/.sort.f90.swp b/prep_regions/.sort.f90.swp
deleted file mode 100644
index eec8c215bcefc974c7ebfaddfc9ff49ee78202f3..0000000000000000000000000000000000000000
Binary files a/prep_regions/.sort.f90.swp and /dev/null differ
diff --git a/prep_syndata/get_error_cov.f90 b/prep_syndata/get_error_cov.f90
index 4131a08da3c99d91a48e67f28fd875b4e4089719..9c05c9a50c9b4f88c81d81795e3196c6be13c167 100644
--- a/prep_syndata/get_error_cov.f90
+++ b/prep_syndata/get_error_cov.f90
@@ -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)
diff --git a/prep_syndata/get_flux.f90 b/prep_syndata/get_flux.f90
index 8ec8d3a667022198c3d173000166cd0e8fdf7864..bbd99509b1009b069fbb1b550516d508bb11ad23 100644
--- a/prep_syndata/get_flux.f90
+++ b/prep_syndata/get_flux.f90
@@ -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
diff --git a/prep_syndata/get_obs.f90 b/prep_syndata/get_obs.f90
index b8661d79a72c684e808fc7d37ee5c943219f625a..3cd38b922e89314a6a17f883377ca84e1193aa0d 100644
--- a/prep_syndata/get_obs.f90
+++ b/prep_syndata/get_obs.f90
@@ -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
diff --git a/prep_syndata/initialize.f90 b/prep_syndata/initialize.f90
index cd6cd6933369ac6d672a083bd10db518a91a919e..42ced70264fb999fc5af182da9d06e0317223395 100644
--- a/prep_syndata/initialize.f90
+++ b/prep_syndata/initialize.f90
@@ -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
diff --git a/prep_syndata/makefile b/prep_syndata/makefile
index c1fc019d9064344a8ed77f083a2413d2618a4d27..827e75ab9bcd6cad951ef676bce2e0b158d2b8a2 100644
--- a/prep_syndata/makefile
+++ b/prep_syndata/makefile
@@ -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)
diff --git a/prep_syndata/mod_perturb.f90 b/prep_syndata/mod_perturb.f90
index e83836e67691c894b4ff220dfee18511a18e56d3..2bdd1b7da16f5e23a4b8aaa70762ca1a84d616ae 100644
--- a/prep_syndata/mod_perturb.f90
+++ b/prep_syndata/mod_perturb.f90
@@ -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
diff --git a/prep_syndata/mod_save.f90 b/prep_syndata/mod_save.f90
index 4b99567791d3fb8766ed9f45e9474f4ab0f6560b..84d34072acf1f2b00a4b8a66f2376aa2f5b7dd72 100644
--- a/prep_syndata/mod_save.f90
+++ b/prep_syndata/mod_save.f90
@@ -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
diff --git a/prep_syndata/mod_settings.f90 b/prep_syndata/mod_settings.f90
index cf2237f79f87b31b1426bdf0379f638ae8034123..479269b226872038988bbf5659450ddbc4570a4d 100644
--- a/prep_syndata/mod_settings.f90
+++ b/prep_syndata/mod_settings.f90
@@ -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
diff --git a/prep_syndata/mod_var.f90 b/prep_syndata/mod_var.f90
index 2aa49892b97975252774d9dfcf926afbaf634eb2..95c68f81949e3d0f536b79983bc4ff74e257158b 100644
--- a/prep_syndata/mod_var.f90
+++ b/prep_syndata/mod_var.f90
@@ -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
! ------------------
diff --git a/prep_syndata/prep_syndata.f90 b/prep_syndata/prep_syndata.f90
index 09e891b9773a2d017129456ce545b505506739f7..ece75880352a873fdfab42f1a5bd05a93b59ef80 100644
--- a/prep_syndata/prep_syndata.f90
+++ b/prep_syndata/prep_syndata.f90
@@ -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
diff --git a/settings/SETTINGS_co2_config b/settings/SETTINGS_co2_config
index 32d38c7c5cd53ee92bedd766669ed6de3d1f956a..2b571e311fdae09d5f6a01a21721da13ff7b3ace 100644
--- a/settings/SETTINGS_co2_config
+++ b/settings/SETTINGS_co2_config
@@ -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
diff --git a/settings/SETTINGS_co2_nest_config b/settings/SETTINGS_co2_nest_config
index 17d270b61aaff10e492c2dbe56a07129e6fe1365..f5237730ea0e6e5fd7cdb9714ff9e9961710985f 100644
--- a/settings/SETTINGS_co2_nest_config
+++ b/settings/SETTINGS_co2_nest_config
@@ -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
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
diff --git a/settings/SETTINGS_ghg_config b/settings/SETTINGS_ghg_config
index 57473da74ab01a2a07887dc3a4bb9fd043b77dc8..05349ab2d0ce59f7e3aba29bc87d80e9df1e5c9d 100644
--- a/settings/SETTINGS_ghg_config
+++ b/settings/SETTINGS_ghg_config
@@ -26,14 +26,17 @@ datef: 20120131
# Inversion method ('analytic', 'congrad' or 'm1qn3')
method: analytic
+# Use lognormal distribution (true or false) (only for non-CO2 species)
+# if true need to use m1qn3 method
+lognormal: .false.
+
+# 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 the optimization (true or false)
inc_ocean: .true.
@@ -58,7 +61,7 @@ restart: 0
# Verbose output
# only use for debugging small runs
-verbose: .true.
+verbose: .false.
# Species ("co2" or "ghg")
spec: ghg
@@ -103,8 +106,8 @@ nt_flx: 12
# used if error in obs input <= zero
measerr: 5.0
-# Initial concentration error: unit same as obs
-cinierr: 5.0
+# Scalar of initial conc error: fraction
+cinierr: 0.01
# Prior flux error: fraction
flxerr: 0.5
diff --git a/settings/SETTINGS_ghg_nest_config b/settings/SETTINGS_ghg_nest_config
index 3d4b28219bf74be06639e8fa9e68f8da57b46315..edbef48446ae44dce8578c7a04cd54d4157c4dde 100644
--- a/settings/SETTINGS_ghg_nest_config
+++ b/settings/SETTINGS_ghg_nest_config
@@ -26,14 +26,17 @@ datef: 20120131
# Inversion method ('analytic', 'congrad' or 'm1qn3')
method: congrad
+# Use lognormal distribution (true or false) (only for non-CO2 species)
+# if true need to use m1qn3 method
+lognormal: .false.
+
+# 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 the optimization (true or false)
inc_ocean: .true.
@@ -103,8 +106,8 @@ nt_flx: 12
# used if error in obs input <= zero
measerr: 5.0
-# Initial concentration error: unit same as obs
-cinierr: 5.0
+# Scalar of initial conc error: fraction
+cinierr: 0.01
# Prior flux error: fraction
flxerr: 0.5
diff --git a/source/error_cov.f90 b/source/error_cov.f90
index 55ac30f89f166a507fba77ceaa2bbe65097b364f..d8a8886340c15db51b254db8677392f704410906 100644
--- a/source/error_cov.f90
+++ b/source/error_cov.f90
@@ -73,6 +73,7 @@ subroutine error_cov(config, files, states, covar, corr, xerr)
real(kind=8), dimension(:), allocatable :: evals
real(kind=8), dimension(:,:), allocatable :: evecs
integer :: ierr, info, neval, ind
+ real, parameter :: smallnum=1.e-15
! calculate error covariance
! --------------------------
@@ -91,14 +92,16 @@ subroutine error_cov(config, files, states, covar, corr, xerr)
end do
covsum = covsum/real(ndt)**2
- ! scale to total error of nested domain (unit Tg/y)
- toterr = sqrt(covsum)*3600.*24.*365./1.e9/numscale
- write(logid,*) 'Total error: ',toterr
- erscalar = config%globerr/toterr
- if ( config%globerr > 0. ) then
- cov = cov * dble(erscalar)**2
- write(logid,*) 'Total error scaled by: ',erscalar
- end if
+ if ( .not.config%lognormal ) then
+ ! scale to total error of nested domain (unit Tg/y)
+ toterr = sqrt(covsum)*3600.*24.*365./1.e9/numscale
+ write(logid,*) 'Total error: ',toterr
+ erscalar = config%globerr/toterr
+ if ( config%globerr > 0. ) then
+ cov = cov * dble(erscalar)**2
+ write(logid,*) 'Total error scaled by: ',erscalar
+ end if
+ endif
! calculate eigendecomposition of cov
! -----------------------------------
@@ -248,7 +251,11 @@ subroutine error_cov(config, files, states, covar, corr, xerr)
! uncertainty for scalars of initial mixing ratio
if ( config%opt_cini ) then
do i = 1, ntcini*ncini
- states%pxerr0(npvar+i) = config%cinierr
+ if ( config%lognormal ) then
+ states%pxerr0(npvar+i) = log(max(smallnum,config%cinierr))
+ else
+ states%pxerr0(npvar+i) = config%cinierr
+ endif
end do
endif
diff --git a/source/get_initcams.f90 b/source/get_initcams.f90
index aff759cbc94d6c314b364ac1b7981b6c4e3186f2..b89a396371a456d102cf43fee0275e7dbdde21ac 100644
--- a/source/get_initcams.f90
+++ b/source/get_initcams.f90
@@ -53,11 +53,13 @@ subroutine get_initcams(filename, varname, ntime, concini)
character(len=*), parameter :: lonname='longitude'
character(len=*), parameter :: latname='latitude'
character(len=*), parameter :: timename='time'
+ character(len=*), parameter :: prsname='hlevel' ! hlevel for CH4, level for N2O
+ logical, parameter :: midpoints=.true. ! true for CH4, false for N2O
integer :: ncid, xid, yid, zid, tid, apid, bpid, varid
- integer :: n, kz
+ integer :: n, kz, status
integer :: nx, ny, nz, nt, nlev, nzmax
integer :: nx_out, ny_out
- logical :: midpoints
+ logical :: flag_hyb
real, dimension(:), allocatable :: lon, lat, lonin, latin
real(kind=8), dimension(:), allocatable :: lono, lato
real, dimension(:), allocatable :: ap, bp, pres, mpres, alt
@@ -66,6 +68,7 @@ subroutine get_initcams(filename, varname, ntime, concini)
real, dimension(nygrid) :: lat_out
real, dimension(:,:,:), allocatable :: work1, work2, work3, work4
real, dimension(:,:,:,:), allocatable :: conc, concin
+ real :: sclfact, offset
! initialize
! ----------
@@ -85,36 +88,60 @@ subroutine get_initcams(filename, varname, ntime, concini)
call check( nf90_open(trim(filename),nf90_NOWRITE,ncid) )
call check( nf90_inq_dimid(ncid,lonname,xid) )
call check( nf90_inq_dimid(ncid,latname,yid) )
- call check( nf90_inq_dimid(ncid,'hlevel',zid) )
+ call check( nf90_inq_dimid(ncid,prsname,zid) )
call check( nf90_inq_dimid(ncid,timename,tid) )
call check( nf90_inquire_dimension(ncid,xid,len=nx) )
call check( nf90_inquire_dimension(ncid,yid,len=ny) )
call check( nf90_inquire_dimension(ncid,zid,len=nlev) )
call check( nf90_inquire_dimension(ncid,tid,len=nt) )
- nz = nlev - 1
if ( ntime.ne.nt ) write(logid,*) 'ERROR get_initcams: inconsistent time dimension'
allocate( lon(nx) )
allocate( lat(ny) )
- allocate( ap(nlev) )
- allocate( bp(nlev) )
- allocate( pres(nlev) )
- allocate( mpres(nz) )
- allocate( alt(nz) )
- allocate( conc(nx,ny,nz,nt) )
! read variables
call check( nf90_inq_varid(ncid,lonname,xid) )
call check( nf90_get_var(ncid,xid,lon) )
call check( nf90_inq_varid(ncid,latname,yid) )
call check( nf90_get_var(ncid,yid,lat) )
- call check( nf90_inq_varid(ncid,'ap',apid) )
- call check( nf90_get_var(ncid,apid,ap) )
- call check( nf90_inq_varid(ncid,'bp',bpid) )
- call check( nf90_get_var(ncid,bpid,bp) )
+ ! check if contains hybrid pressure coordinates
+ status = nf90_inq_varid(ncid,'hyai',apid)
+ print*, 'get_initcams: status = ',status
+ if (status.eq.-49) then
+ ! no hybrid pressure coordinates
+ flag_hyb = .false.
+ nz = nlev
+ allocate( mpres(nz) )
+ call check( nf90_inq_varid(ncid,prsname,zid) )
+ call check( nf90_get_var(ncid,zid,mpres) )
+ else
+ ! contains hybrid pressure coordinates
+ flag_hyb = .true.
+ nz = nlev - 1
+ allocate( ap(nlev) )
+ allocate( bp(nlev) )
+ allocate( pres(nlev) )
+ allocate( mpres(nz) )
+ call check( nf90_inq_varid(ncid,'hyai',apid) )
+ call check( nf90_get_var(ncid,apid,ap) )
+ call check( nf90_inq_varid(ncid,'hybi',bpid) )
+ call check( nf90_get_var(ncid,bpid,bp) )
+ endif
+ allocate( alt(nz) )
+ allocate( conc(nx,ny,nz,nt) )
call check( nf90_inq_varid(ncid,trim(varname),varid) )
call check( nf90_get_var(ncid,varid,conc(:,:,:,:)) )
+ ! check if file contains scalar and offset
+ status = nf90_inquire_attribute(ncid,varid,'scale_factor')
+ if ( status.eq.nf90_noerr ) then
+ status = nf90_inquire_attribute(ncid,varid,'add_offset')
+ if ( status.eq.nf90_noerr ) then
+ call check( nf90_get_att(ncid, varid, 'scale_factor', sclfact) )
+ call check( nf90_get_att(ncid, varid, 'add_offset', offset) )
+ conc = conc*sclfact + offset
+ endif
+ endif
call check( nf90_close(ncid) )
! adjust lat and lon of grid
@@ -124,9 +151,13 @@ subroutine get_initcams(filename, varname, ntime, concini)
allocate( lono(nx) )
lato = dble(lat)
lono = dble(lon)
+ ! sort from south to north and west to east
call sort(ny, lato, indy)
call sort(nx, lono, indx)
conc = conc(indx,indy,:,:)
+ print*, 'get_initcams: indy = ',indy
+ print*, 'get_initcams: indx = ',indx
+ ! cut extra longitude grid cells (assumes grid is -180 to 180 and not 0 to 360)
if ( lono(nx).ge.180. ) then
allocate( concin((nx-1),ny,nz,nt) )
allocate( lonin((nx-1)) )
@@ -141,26 +172,36 @@ subroutine get_initcams(filename, varname, ntime, concini)
allocate( latin(ny) )
if ( lato(ny).eq.90. ) then
do n = 1, ny
- latin(n) = -90. + real(n - 1)*(180./96.)
+ latin(n) = -90. + real(n - 1)*(180./real(ny))
end do
else
latin = sngl(lato)
endif
-! print*, 'get_initcams: latin = ',latin
-! print*, 'get_initcams: lonin = ',lonin
+ print*, 'get_initcams: latin = ',latin
+ print*, 'get_initcams: lonin = ',lonin
+
+ ! calculate mid layer pressure
+ ! (for simplicity currently use standard surface pressure
+ ! should be good enough for background calculation,
+ ! future update would be to read surf pres from file and do
+ ! vertical interpolation to flexpart levels grid by grid)
+ if ( flag_hyb ) then
+ do n = 1, nlev
+ pres(n) = ap(n) + bp(n)*psurf
+ end do
+ ! avoid pressure of zero
+ pres(nlev) = 2.
+ do n = 1, nz
+ mpres(n) = 0.5*(pres(n) + pres(n+1))
+ end do
+ endif
- ! calculate altitude
- do n = 1, nlev
- pres(n) = ap(n) + bp(n)*psurf
- end do
- ! avoid pressure of zero
- pres(nlev) = 2.
do n = 1, nz
- mpres(n) = 0.5*(pres(n) + pres(n+1))
- alt(n)=7500*log(psurf/pres(n+1))
+ alt(n)=7500*log(psurf/mpres(n))
end do
-! print*, 'get_initcams: pres = ',pres
-! print*, 'get_initcams: alt = ',alt
+ if ( flag_hyb ) print*, 'get_initcams: pres = ',pres
+ print*, 'get_initcams: mpres = ',mpres
+ print*, 'get_initcams: alt = ',alt
! interpolate to flexpart grid
! ----------------------------
@@ -175,13 +216,12 @@ subroutine get_initcams(filename, varname, ntime, concini)
! determine upper alt level of conc field for vertical interpolation
nzmax = minloc(alt, dim=1, mask=abs(alt-outheight(nzgrid)).eq.minval(abs(alt-outheight(nzgrid))))
nzmax = min(nz, nzmax+1)
-! print*, 'get_initcams: nzmax = ',nzmax
+ print*, 'get_initcams: nzmax = ',nzmax
allocate( work2(nxgrid,nygrid,1) )
allocate( work3(nxgrid,nygrid,nzmax) )
allocate( work4(nxgrid,nygrid,nzgrid) )
- midpoints = .false.
do n = 1, nt
! interpolate horizontally
do kz = 1, nzmax
@@ -203,9 +243,9 @@ subroutine get_initcams(filename, varname, ntime, concini)
deallocate( latin )
deallocate( indx )
deallocate( indy )
- deallocate( ap )
- deallocate( bp )
- deallocate( pres )
+ if ( allocated(ap) ) deallocate( ap )
+ if ( allocated(bp) ) deallocate( bp )
+ if ( allocated(pres) ) deallocate( pres )
deallocate( mpres )
deallocate( alt )
deallocate( conc )
diff --git a/source/init_cams.f90 b/source/init_cams.f90
index 0f9121ff81c0f9c06688e669712bdd14ca2b53d7..4d6f33a952a8363ee79ce4b54136f904c0d08b52 100644
--- a/source/init_cams.f90
+++ b/source/init_cams.f90
@@ -56,8 +56,11 @@ subroutine init_cams(files, config, obs)
character(max_path_len) :: path_flexrec, filename
character(max_name_len) :: varname
character(len=8) :: adate, areldate, areltime
+ character(len=4) :: ayear
+ character(len=2) :: amonth
logical :: lexist
+ integer, parameter :: hfreq = 24 ! CAMS N2O is 3-hourly but CH4 is 24-hourly
integer :: jjjjmm, jjjjmmdd, hhmiss, eomday
integer :: lasttime, startmonth
integer :: month, prevmonth
@@ -99,6 +102,9 @@ subroutine init_cams(files, config, obs)
month = jjjjmmdd/100
path_flexrec = trim(files%path_flexpart)//trim(obs%recs(ind(i)))//'/'//trim(adate)//'/'
+ print*, 'init_cams: path_flexrec = ',path_flexrec
+ print*, 'init_cams: average_fp = ',config%average_fp
+
if ( config%average_fp ) then
! match and, if needed, average grid_initial files to observation timestamp
@@ -187,15 +193,15 @@ subroutine init_cams(files, config, obs)
call caldate(jdi, jjjjmmdd, hhmiss)
jjjjmm = jjjjmmdd/100
! round time to nearest 3-hours
- hh = (jdi - floor(jdi))*24d0
- hh = floor(hh/3d0)*3d0
-
+! hh = (jdi - floor(jdi))*24d0
+! hh = floor(hh/3d0)*3d0
+
! read concentration field only if not already in memory
if ( jjjjmm.ne.lasttime ) then
- ! allocate concini knowing data are 3-hourly for each month
+ ! allocate concini knowing data freq (from hfreq)
if (allocated( concini ) ) deallocate( concini )
eomday = calceomday(jjjjmm)
- ntime = eomday*8
+ ntime = eomday*24/hfreq
allocate( concini(nxgrid,nygrid,nzgrid,ntime), stat = ierr )
if ( ierr.ne.0 ) then
write(logid,*) 'ERROR init_cams: not enough memory'
@@ -203,9 +209,10 @@ subroutine init_cams(files, config, obs)
endif
concini(:,:,:,:) = 0.
! read new concentration field
- write(adate,'(I6.6)') jjjjmm
- print*, 'init_cams: file initial = ',files%file_initconc
- filename = str_replace(files%file_initconc, 'YYYYMM', adate)
+ write(ayear,'(I4.4)') jjjjmm/100
+ write(amonth,'(I2.2)') jjjjmm-(jjjjmm/100)*100
+ filename = str_replace(files%file_initconc, 'YYYY', ayear)
+ filename = str_replace(filename, 'MM', amonth)
filename = trim(files%path_initconc)//trim(filename)
print*, 'init_cams: filename = ',filename
inquire(file=trim(filename),exist=lexist)
@@ -213,27 +220,17 @@ subroutine init_cams(files, config, obs)
varname = files%varname_init
call get_initcams(filename, varname, ntime, concini)
else
- ! if doesn't exist try using first file of the year
- write(adate,'(I6.6)') (jjjjmm/100)*100+1
- filename = str_replace(files%file_initconc, 'YYYYMM', adate)
- filename = trim(files%path_initconc)//trim(filename)
- print*, 'init_cams: filename = ',filename
- inquire(file=trim(filename),exist=lexist)
- if(lexist) then
- varname = files%varname_init
- call get_initcams(filename, varname, ntime, concini)
- else
- write(logid,*) 'WARNING: cannot find ',filename
- go to 10
- endif
+ write(logid,*) 'ERROR init_cams: cannot find file ',filename
+ stop
endif
endif
! calculate contribution of initial concentration to receptor
! index time dimension of concini
+ ! for CAMS N2O 3-hourly data, for CAMS CH4 daily data
startmonth = (jjjjmmdd/100)*100 + 1
nhours = int((jdi - juldate(startmonth, 0))*24d0)
- nt = nhours/3 + 1
+ nt = nhours/hfreq + 1
print*, 'init_cams: startmonth = ',startmonth
print*, 'init_cams: nhours = ',nhours
print*, 'init_cams: nt = ',nt
diff --git a/source/init_ghg.f90 b/source/init_ghg.f90
index b4a9fd7da67332845fe5054b4ae1fcd186ae4631..e6f233a4682f1c2e9a644d56aede893eabb1590b 100644
--- a/source/init_ghg.f90
+++ b/source/init_ghg.f90
@@ -90,6 +90,7 @@ subroutine init_ghg(files, config, fluxes, obs, states, covar)
character(len=3), dimension(:), allocatable:: recs
real :: conc
character(len=20) :: dump
+ real, parameter :: smallnum=1.e-15
! global prior fluxes
! -------------------
@@ -112,7 +113,7 @@ subroutine init_ghg(files, config, fluxes, obs, states, covar)
call caldate(jd, yyyymmdd, hhmiss)
write(adate,'(I4.4)') yyyymmdd/10000
- nread = int(statres/(real(365)/real(config%nt_flx)))
+ nread = nint(statres/(real(365)/real(config%nt_flx)))
if ( nread.eq.0 ) nread = 1
allocate( flx(nxgrid,nygrid,nread), stat = ierr )
if ( ierr.ne.0 ) stop 'ERROR init_ghg: not enough memory'
@@ -224,16 +225,10 @@ subroutine init_ghg(files, config, fluxes, obs, states, covar)
do ix = 1, nxregrid
do jy = 1, nyregrid
area = grid_area(reg_lat(jy)+0.5*rdy, rdy, rdx)
+ ! if inc_ocean flx_box contains all fluxes, if false only land fluxes
if ( nbox_xy(ix,jy).gt.0 ) then
flx_box(nbox_xy(ix,jy),n) = flx_box(nbox_xy(ix,jy),n) + &
flx(ix0+ix-1,jy0+jy-1,n) * area
- endif
- if ( (nbox_xy(ix,jy).gt.0).and.(lsm(ix,jy).ge.0.01) ) then
- ! land box -> error is proportional to flux
- err_box(nbox_xy(ix,jy),n) = err_box(nbox_xy(ix,jy),n) + &
- flx(ix0+ix-1,jy0+jy-1,n) * config%flxerr * area
- else if ( config%inc_ocean ) then
- ! ocean box -> error is proportional to flux
err_box(nbox_xy(ix,jy),n) = err_box(nbox_xy(ix,jy),n) + &
flx(ix0+ix-1,jy0+jy-1,n) * config%flxerr * area
endif
@@ -289,11 +284,36 @@ subroutine init_ghg(files, config, fluxes, obs, states, covar)
! allocate state and error vectors
call alloc_states(states)
- ! prior state vector zero since optimize offsets
- states%px0(1:npvar) = 0.
+ if ( config%lognormal ) then
+ ! lognormal distribution
+! ! replace values <= 0 with smallnum since need to transform to log
+! do n = 1, ntstate
+! where( flx_box(:,n).lt.smallnum ) flx_box(:,n) = smallnum
+! states%px0((n-1)*nbox+1:n*nbox) = flx_box(:,n)
+! end do
+ ! optimize scalar z = ln(x/xb)
+ ! prior state vector z is zero
+ states%px0(1:npvar) = 0.
+ else
+ ! normal distribution
+ ! prior state vector zero since optimize offsets
+ states%px0(1:npvar) = 0.
+ endif
! assign prior initial mixing ratio scalars
- if ( config%opt_cini ) states%px0(npvar+1:nvar) = 1.
+ if ( config%opt_cini ) then
+ if ( config%lognormal ) then
+! states%px0(npvar+1:nvar) = 1.
+ states%px0(npvar+1:nvar) = 0.
+ else
+ states%px0(npvar+1:nvar) = 1.
+ endif
+ endif
+
+ ! lognormal transform to log
+! if ( config%lognormal ) then
+! states%px0(:) = log(states%px0)
+! endif
! time stamp of flux variables
states%xtime(:) = fluxes%time
@@ -304,8 +324,19 @@ subroutine init_ghg(files, config, fluxes, obs, states, covar)
end do
print*, 'init_ghg: cinitime = ',states%cinitime
- ! average error over time steps
- xerr = abs(sum(err_box(:,:),dim=2)/real(ntstate) )
+ if ( config%lognormal ) then
+ ! lognormal
+! ! errors are for log transformed state variables
+! xerr = sum(flx_box(:,:),dim=2)/real(ntstate)
+! xerr = log(xerr)
+! xerr = config%flxerr*xerr
+ ! uncertainty of a scalar of fluxes -> use constant uncertainty
+ xerr(:) = config%flxerr
+ else
+ ! normal
+ ! average error over time steps
+ xerr = abs(sum(err_box(:,:),dim=2)/real(ntstate) )
+ endif
! assign best guess of state vector
if ( (config%prior_bg.eq.1).and.(config%method.eq.'congrad') ) then
@@ -316,6 +347,13 @@ subroutine init_ghg(files, config, fluxes, obs, states, covar)
states%px(:) = states%px0(:)
endif
+ ! testing
+ open(100,file=trim(files%path_output)//'prior.txt',status='replace',action='write',iostat=ierr)
+ do n = 1, nvar
+ write(100,fmt='(E14.8)') states%px(n)
+ end do
+ close(100)
+
! correlation and error covariance matrices
! -----------------------------------------
@@ -355,6 +393,7 @@ subroutine init_ghg(files, config, fluxes, obs, states, covar)
! initial concentrations from termination of trajectories
! only do for a fresh run
if ( config%restart.eq.0 ) then
+! call init_cams(files, config, obs)
call init_cini(files, config, obs)
! call init_fpctm(files, config, obs)
endif
diff --git a/source/m1qn3_interface.f90 b/source/m1qn3_interface.f90
index e81417189122e2a83d239f7bd4acfda4710576e7..52a3d1927ab804e4379f26bd453a2b1c65cae1d9 100644
--- a/source/m1qn3_interface.f90
+++ b/source/m1qn3_interface.f90
@@ -67,7 +67,7 @@ subroutine m1qn3_interface(iter, grad, files, config, fluxes, obs, states, covar
! local variables for interface with simulate.f90
real :: cost_o, cost_p
- real(kind=8), dimension(nvar,1) :: phi, chi_ad, px
+ real(kind=8), dimension(nvar,1) :: phi, chi_ad, px, ones
real, dimension(nvar) :: grad_o
character(len=2) :: aiter
@@ -89,6 +89,7 @@ subroutine m1qn3_interface(iter, grad, files, config, fluxes, obs, states, covar
! initialize local variables
! --------------------------
+ ones(:,1) = 1d0
n = nvar
! initial guess at vector that minimizes the cost in chi space
x(:) = 0d0
@@ -186,8 +187,15 @@ subroutine m1qn3_interface(iter, grad, files, config, fluxes, obs, states, covar
grad_o(:) = 0.
call simulate(iter, files, config, fluxes, obs, states, grad_o, cost_o)
! calculate cost
- ! Jp = x^Tx where x = state vector in chi space
- cost_p = dot_product(x,x)
+ if ( config%lognormal ) then
+ ! lognormal
+ ! 2Jp = 2*sum(z) + x^Tx (where x = state vector in chi space and z is state vector in physical space)
+ cost_p = sum(states%px) + dot_product(x,x)
+ else
+ ! normal
+ ! 2Jp = x^Tx (where x = state vector in chi space)
+ cost_p = dot_product(x,x)
+ endif
! J = 0.5*(Jp + Jo)
cost = 0.5*(cost_p + cost_o)
f = cost
@@ -195,7 +203,7 @@ subroutine m1qn3_interface(iter, grad, files, config, fluxes, obs, states, covar
write(logid,*) 'Cost iteration '//aiter//': ',cost
! save cost in case of a restart
open(100,file=trim(files%path_output)//'cost.txt',status='unknown',access='append',action='write',iostat=ierr)
- write(100,fmt='(F11.2)') cost
+ write(100,fmt='(F13.4)') cost
close(100)
endif
@@ -203,8 +211,18 @@ subroutine m1qn3_interface(iter, grad, files, config, fluxes, obs, states, covar
phi(:,1) = dble(grad_o)
chi_ad = phi2chi_ad(config, covar, phi )
- ! new gradient J' = J'p + J'o
- g(:) = x(:) + chi_ad(:,1)
+ ! full gradient chi space
+ ! J' = J'p + J'o
+ if ( config%lognormal ) then
+ ! lognormal
+ ! J'p = B^(1/2)ones + chi
+ ! misuse phi2chi_ad to calculate B^(1/2)ones
+ phi = phi2chi_ad(config, covar, ones )
+ g(:) = phi(:,1) + x(:) + chi_ad(:,1)
+ else
+ ! normal
+ g(:) = x(:) + chi_ad(:,1)
+ endif
end do iterate
diff --git a/source/makefile b/source/makefile
index 7c777b157396024085a249370fcbb6937c608f58..103ee50dabf89f195ce7a7e982c33907a740b179 100644
--- a/source/makefile
+++ b/source/makefile
@@ -4,12 +4,12 @@ INCPATH = /usr/include/
LNK = -o
CMPL = -c
LIBS = -lnetcdf -lnetcdff -llapack
-FFLAGS = -O2 -m64 -mcmodel=medium -fconvert=little-endian -frecord-marker=4 -ffree-form -ffree-line-length-0
-#FFLAGS = -O0 -g -m64 -fbounds-check -mcmodel=medium -fconvert=little-endian -frecord-marker=4 -ffree-form \
-# -fbacktrace -ffree-line-length-0 -fcheck=all -Wall
+#FFLAGS = -O2 -m64 -mcmodel=medium -fconvert=little-endian -frecord-marker=4 -ffree-form -ffree-line-length-0
+FFLAGS = -O0 -g -m64 -fbounds-check -mcmodel=medium -fconvert=little-endian -frecord-marker=4 -ffree-form \
+ -fbacktrace -ffree-line-length-0 -fcheck=all -Wall
LDFLAGS = $(FFLAGS) -L$(LIBPATH) -I$(INCPATH) $(LIBS)
-MAIN = main_test
+MAIN = main
SRCS = mod_var.f90 \
mod_settings.f90 \
diff --git a/source/mod_perturb.f90 b/source/mod_perturb.f90
index 133b6f6257a09ae6e8695f7884e77193261bfb54..23beee384983809fc69fc0715e428bb8b1cd013d 100644
--- a/source/mod_perturb.f90
+++ b/source/mod_perturb.f90
@@ -96,7 +96,14 @@ module mod_perturb
endif
! add perturbation to state vector
- states%px0(:) = states%px0(:) + real(perturb(:),kind=4)
+ if ( config%lognormal ) then
+ ! lognormal
+! states%px0(:) = states%px0(:) * exp(real(perturb(:),kind=4))
+ states%px0(:) = exp(real(perturb(:),kind=4))
+ else
+ ! normal
+ states%px0(:) = states%px0(:) + real(perturb(:),kind=4)
+ endif
if ( config%prior_bg.ne.1 ) then
states%px(:) = states%px0(:)
endif
diff --git a/source/mod_save.f90 b/source/mod_save.f90
index 571f49627d2a76946f7c6bd37efc5c61affbedc4..8d139246f636afa1b8d2b18de8a3e093183401b0 100644
--- a/source/mod_save.f90
+++ b/source/mod_save.f90
@@ -221,7 +221,7 @@ module mod_save
end do
end do
- ! redistribute so that flux increments on ocean grid cells are attributed to
+ ! redistribute so that flux increments on ocean boxes are attributed to
! land grid cells within the same aggregated box
! note: future improvement to define aggregated grid with clear separation along coastlines
! then this step no longer needed and should be more accurate!
@@ -451,30 +451,85 @@ module mod_save
do n = 1, ntstate
do jy = 1, nyregrid
do ix = 1, nxregrid
+ if ( config%nested ) then
+ fpri(ix,jy,n) = fluxes%flx_nest(ix,jy,n)/numscale
+ else
+ fpri(ix,jy,n) = fluxes%flx(ix0+ix-1,jy0+jy-1,n)/numscale
+ endif
+ ! 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 ( config%nested ) then
- fpri(ix,jy,n) = fluxes%flx_nest(ix,jy,n)/numscale
- else
- fpri(ix,jy,n) = fluxes%flx(ix0+ix-1,jy0+jy-1,n)/numscale
- endif
! regrid flux increments according to prior flux distribution and add prior flux
if ( fluxes%flx_box(nbox_xy(ix,jy),n,1).gt.smallnum ) then
if ( config%nested ) then
- fpos(ix,jy,n) = states%px((n-1)*nbox+nbox_xy(ix,jy)) * &
- fluxes%flx_nest(ix,jy,n)/fluxes%flx_box(nbox_xy(ix,jy),n,1)/numscale + &
- fpri(ix,jy,n)
+ if ( config%lognormal ) then
+! ! lognormal distribution: optimize fluxes not offsets
+! fpos(ix,jy,n) = exp(states%px((n-1)*nbox+nbox_xy(ix,jy))) * &
+! fluxes%flx_nest(ix,jy,n)/fluxes%flx_box(nbox_xy(ix,jy),n,1)/numscale
+! ! optimize scalars of fluxes
+ fpos(ix,jy,n) = exp(states%px((n-1)*nbox+nbox_xy(ix,jy))) * fluxes%flx_nest(ix,jy,n)/numscale
+ else
+ ! normal distribution: optimize offsets
+ fpos(ix,jy,n) = states%px((n-1)*nbox+nbox_xy(ix,jy)) * &
+ fluxes%flx_nest(ix,jy,n)/fluxes%flx_box(nbox_xy(ix,jy),n,1)/numscale + &
+ fpri(ix,jy,n)
+ endif
else
- fpos(ix,jy,n) = states%px((n-1)*nbox+nbox_xy(ix,jy)) * &
- fluxes%flx(ix0+ix-1,jy0+jy-1,n)/fluxes%flx_box(nbox_xy(ix,jy),n,1)/numscale + &
- fpri(ix,jy,n)
+ if ( config%lognormal ) then
+ ! lognormal distribution: optimize fluxes not offsets
+! fpos(ix,jy,n) = exp(states%px((n-1)*nbox+nbox_xy(ix,jy))) * &
+! fluxes%flx(ix0+ix-1,jy0+jy-1,n)/fluxes%flx_box(nbox_xy(ix,jy),n,1)/numscale
+! ! optimize scalars of fluxes
+ fpos(ix,jy,n) = exp(states%px((n-1)*nbox+nbox_xy(ix,jy))) * fluxes%flx(ix0+ix-1,jy0+jy-1,n)/numscale
+ else
+ ! normal distribution: optimize offsets
+ fpos(ix,jy,n) = states%px((n-1)*nbox+nbox_xy(ix,jy)) * &
+ fluxes%flx(ix0+ix-1,jy0+jy-1,n)/fluxes%flx_box(nbox_xy(ix,jy),n,1)/numscale + &
+ fpri(ix,jy,n)
+ endif
endif
else
- fpos(ix,jy,n) = states%px((n-1)*nbox+nbox_xy(ix,jy))/numscale + fpri(ix,jy,n)
+ if ( config%lognormal ) then
+! ! lognormal distribution: optimize fluxes not offsets
+! fpos(ix,jy,n) = exp(states%px((n-1)*nbox+nbox_xy(ix,jy)))/numscale
+! ! optimize scalars of fluxes
+ if ( config%nested ) then
+ fpos(ix,jy,n) = exp(states%px((n-1)*nbox+nbox_xy(ix,jy))) * fluxes%flx_nest(ix,jy,n)/numscale
+ else
+ fpos(ix,jy,n) = exp(states%px((n-1)*nbox+nbox_xy(ix,jy))) * fluxes%flx(ix0+ix-1,jy0+jy-1,n)/numscale
+ endif
+ else
+ ! normal distribution: optimize offsets
+ fpos(ix,jy,n) = states%px((n-1)*nbox+nbox_xy(ix,jy))/numscale + fpri(ix,jy,n)
+ endif
endif
! regrid without adding prior fluxes
- epri(ix,jy,n) = states%pxerr0((n-1)*nbox+nbox_xy(ix,jy))/numscale
- epos(ix,jy,n) = states%pxerr((n-1)*nbox+nbox_xy(ix,jy))/numscale
- xpos(ix,jy,n) = states%px((n-1)*nbox+nbox_xy(ix,jy))/numscale
+ if ( config%lognormal ) then
+ ! lognormal distribution
+! epri(ix,jy,n) = exp(states%pxerr0((n-1)*nbox+nbox_xy(ix,jy)))/numscale
+! epos(ix,jy,n) = exp(states%pxerr((n-1)*nbox+nbox_xy(ix,jy)))/numscale
+! xpos(ix,jy,n) = exp(states%px((n-1)*nbox+nbox_xy(ix,jy)))/numscale
+ ! errors are for z = ln(x/xb) so errors in x = xb*(exp(err) - 1)
+ if ( config%nested ) then
+ epri(ix,jy,n) = fluxes%flx_nest(ix,jy,n)*abs(exp(states%pxerr0((n-1)*nbox+nbox_xy(ix,jy))) - 1.)/numscale
+ epos(ix,jy,n) = fluxes%flx_nest(ix,jy,n)*abs(exp(states%pxerr((n-1)*nbox+nbox_xy(ix,jy))) - 1.)/numscale
+ else
+ epri(ix,jy,n) = fluxes%flx(ix0+ix-1,jy0+jy-1,n)*abs(exp(states%pxerr0((n-1)*nbox+nbox_xy(ix,jy))) - 1.)/numscale
+ epos(ix,jy,n) = fluxes%flx(ix0+ix-1,jy0+jy-1,n)*abs(exp(states%pxerr((n-1)*nbox+nbox_xy(ix,jy))) - 1.)/numscale
+ endif
+ xpos(ix,jy,n) = exp(states%px((n-1)*nbox+nbox_xy(ix,jy)))
+ else
+ ! normal distribution
+ epri(ix,jy,n) = states%pxerr0((n-1)*nbox+nbox_xy(ix,jy))/numscale
+ epos(ix,jy,n) = states%pxerr((n-1)*nbox+nbox_xy(ix,jy))/numscale
+ xpos(ix,jy,n) = states%px((n-1)*nbox+nbox_xy(ix,jy))/numscale
+ endif
+ else
+ ! add ocean fluxes if not optimized
+ if ( config%nested ) then
+ fpos(ix,jy,n) = fluxes%flx_nest(ix,jy,n)/numscale
+ else
+ fpos(ix,jy,n) = fluxes%flx(ix0+ix-1,jy0+jy-1,n)/numscale
+ endif
endif
end do
end do
@@ -484,19 +539,27 @@ module mod_save
totpos = 0.
do n = 1, ntstate
do nb = 1, nbox
- totpos = totpos + states%px((n-1)*nbox+nb)*area_box(nb)/numscale
+ if ( config%lognormal ) then
+ totpos = totpos + exp(states%px((n-1)*nbox+nb))*fluxes%flx_box(nb,n,1)*area_box(nb)/numscale
+ else
+ totpos = totpos + states%px((n-1)*nbox+nb)*area_box(nb)/numscale
+ endif
end do
end do
totpos = totpos/real(ntstate)
- print*, 'Post total before regridding = ',totpos*3600.*24.*365./1.e9
+ print*, 'Total increment before regridding = ',totpos*3600.*24.*365./1.e9
totpos = 0.
do jy = 1, nyregrid
area = grid_area(reg_lat(jy)+0.5*rdy, rdy, rdx)
do ix = 1, nxregrid
- totpos = totpos + sum(fpos(ix,jy,:)-fpri(ix,jy,:))*area/real(ntstate)
+ if ( config%lognormal ) then
+ totpos = totpos + sum(fpos(ix,jy,:))*area/real(ntstate)
+ else
+ totpos = totpos + sum(fpos(ix,jy,:)-fpri(ix,jy,:))*area/real(ntstate)
+ endif
end do
end do
- print*, 'Post total after regridding = ',totpos*3600.*24.*365./1.e9
+ print*, 'Total increment after regridding = ',totpos*3600.*24.*365./1.e9
! write ncdf file
! ---------------
diff --git a/source/mod_settings.f90 b/source/mod_settings.f90
index baaeeeda61f2ed9d78ebc28824e0dc15c9f7cc77..a7942ecd1404f93ed5364f5d9564617d34232b44 100644
--- a/source/mod_settings.f90
+++ b/source/mod_settings.f90
@@ -109,7 +109,8 @@ module mod_settings
integer :: datef ! end date (yyyymmdd)
logical :: average_fp ! average footprints to match observation (true or false)
character(len=3) :: spec ! species ('co2' or 'ghg')
- character(len=max_name_len) :: method ! inversion method ('analytic' or 'congrad')
+ character(len=max_name_len) :: method ! inversion method ('analytic' or 'congrad' or 'm1qn3')
+ logical :: lognormal ! use log-normal distribution in state space
integer :: maxiter ! max number of iterations for congrad and m1qn3
logical :: inc_ocean ! include ocean boxes in inversion (true or false)
logical :: opt_cini ! optimize background mixing ratios (true or false)
@@ -607,6 +608,7 @@ module mod_settings
config%opt_cini = .false.
config%spa_corr = .false.
config%verbose = .false.
+ config%lognormal = .false.
! open file
open (100, file = trim (filename), status = 'old', iostat=ierr)
@@ -656,6 +658,9 @@ module mod_settings
identifier = "method:"
call read_content (line, identifier, cc, cn, cl, match)
if ( match ) config%method = cc
+ identifier = "lognormal:"
+ call read_content (line, identifier, cc, cn, cl, match)
+ if ( match ) config%lognormal = cl
identifier = "maxiter:"
call read_content (line, identifier, cc, cn, cl, match)
if ( match ) config%maxiter = int(cn)
@@ -776,41 +781,53 @@ module mod_settings
! 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'
+ print*, 'WARNING read_config_settings: background optimization should not be used for CO2 -> switching off'
config%opt_cini = .false.
endif
! initial checks
if ( config%spa_corr ) then
if ( config%sigma_land.lt.0 ) then
- write(logid,*) 'ERROR: uncorrect value of sigma_land'
+ print*, 'ERROR: uncorrect value of sigma_land'
stop
endif
if ( config%sigma_ocean.lt.0 ) then
- write(logid,*) 'ERROR: uncorrect value of sigma_ocean'
+ print*, 'ERROR: uncorrect value of sigma_ocean'
stop
endif
if ( config%sigmatime.lt.0 ) then
- write(logid,*) 'ERROR: uncorrect value of sigmatime'
+ print*, 'ERROR: uncorrect value of sigmatime'
stop
endif
endif
if ( config%measerr.lt.0 ) then
- write(logid,*) 'ERROR: uncorrect value of measerr'
+ print*, 'ERROR: uncorrect value of measerr'
stop
endif
if ( config%cinierr.lt.0 ) then
- write(logid,*) 'ERROR: uncorrect value of cinierr'
+ print*, 'ERROR: uncorrect value of cinierr'
stop
endif
if ( config%flxerr.lt.0 ) then
- write(logid,*) 'ERROR: uncorrect value of flxerr'
+ print*, 'ERROR: uncorrect value of flxerr'
stop
endif
if ( config%flxerr_ll.lt.0 ) then
- write(logid,*) 'ERROR: uncorrect value of flxerr_ll'
+ print*, 'ERROR: uncorrect value of flxerr_ll'
stop
endif
+ if ( config%lognormal ) then
+ print*, 'MESSAGE: lognormal inversion needs M1QN3 method'
+ config%method = 'm1qn3'
+ if ( config%spec.eq.'co2' ) then
+ print*, 'WARNING read_config_settings: for CO2 cannot use lognormal -> changing to normal'
+ config%lognormal = .false.
+ endif
+ endif
+ if ( config%spec.eq.'co2'.and.config%inc_ocean ) then
+ print*, 'WARNING read_config_settings: for CO2 ocean fluxes are not optimized inc_ocean -> false'
+ config%inc_ocean = .false.
+ endif
end subroutine read_config_settings
diff --git a/source/read_obs.f90 b/source/read_obs.f90
index 60e054ef0c8e87918734e018eb5394c628905b8d..e1b276071f5af64189af62408f0f946520acf863 100755
--- a/source/read_obs.f90
+++ b/source/read_obs.f90
@@ -123,7 +123,7 @@ subroutine read_obs(config, files)
read(100,fmt='(A3,1X,I8,1X,I6,1X,F14.6,1X,F10.6,1X,F10.4,1X,F10.4,1X,I4)',iostat=ierr) &
recs, yyyymmdd, hhmiss, jdate, avetime, conc, err, num
endif
- ! if avetime not included (backwards compatability with observations prepared before update)
+ ! if avetime not included (backwards compatability with observations prepared before update)
! avetime = 0d0
! if (reclen.eq.4) then
! read(100,fmt='(A4,1X,I8,1X,I6,1X,F14.6,1X,F10.4,1X,F10.4,1X,I4)',iostat=ierr) &
diff --git a/source/retrieval.f90 b/source/retrieval.f90
index ffe305198ae7be81bdbd39c529bd5a548d58232e..93be9898226634d6442a5fe7b5c499ee2b6fda83 100644
--- a/source/retrieval.f90
+++ b/source/retrieval.f90
@@ -40,6 +40,7 @@
!! congrad
!! analytic
!! prod_icov
+!! m1qn3_interface
!!
!---------------------------------------------------------------------------------------
@@ -63,8 +64,8 @@ subroutine retrieval(files, config, fluxes, obs, states, covar)
type (states_t), intent (in out) :: states
type (covar_t), intent (in out) :: covar
- real, dimension(nvar) :: grad_p, grad_o
- real(kind=8), dimension(nvar,1) :: grad, chi_ad
+ real, dimension(nvar) :: grad_p, grad_o, icovp
+ real(kind=8), dimension(nvar,1) :: grad, chi_ad, phi, ones, chi, chi1
real :: cost_p, cost_o, cost
integer :: iter
integer :: ierr, n
@@ -78,6 +79,7 @@ subroutine retrieval(files, config, fluxes, obs, states, covar)
grad(:,1) = 0d0
cost_o = 0.
cost_p = 0.
+ ones(:,1) = 1d0
if ( config%method.eq.'analytic' ) then
! allocate transport matrix
@@ -97,12 +99,27 @@ subroutine retrieval(files, config, fluxes, obs, states, covar)
! -------------------------
! gradient state space
- ! J'p = B^(-1)(p - p0)
- grad_p = prod_icov(config, covar, states%px - states%px0)
+ if ( config%lognormal ) then
+ ! lognormal (note p is log transformed variable)
+ ! J'p = ones + B^(-1)(p - p0)
+ grad_p = ones(:,1) + prod_icov(config, covar, states%px - states%px0)
+ else
+ ! normal
+ ! J'p = B^(-1)(p - p0)
+ grad_p = prod_icov(config, covar, states%px - states%px0)
+ endif
! cost state space
- ! 2Jp = (p - p0)^TB^(-1)(p - p0)
- cost_p = dot_product( states%px - states%px0, real(grad_p,kind=4) )
+ if ( config%lognormal ) then
+ ! lognormal (note p is log transformed variable)
+ ! 2Jp = 2*sum(p) + (p - p0)^TB^(-1)(p - p0)
+ icovp = prod_icov(config, covar, states%px - states%px0)
+ cost_p = 2.*sum(states%px) + dot_product( states%px - states%px0, icovp )
+ else
+ ! normal
+ ! 2Jp = (p - p0)^TB^(-1)(p - p0)
+ cost_p = dot_product( states%px - states%px0, grad_p )
+ endif
! gradient observation space
! J'o = H^(T)R^(-1)(H(p) - y)
@@ -154,7 +171,7 @@ subroutine retrieval(files, config, fluxes, obs, states, covar)
call simulate(iter, files, config, fluxes, obs, states, grad_o, cost_o)
else if ( config%method.eq.'m1qn3' ) then
! total gradient
- ! J' = J'p + J'o
+ ! J' = J'p + J'o
grad(:,1) = dble(grad_p + grad_o)
! transform grad to chi space
chi_ad = phi2chi_ad(config, covar, grad)
@@ -169,9 +186,17 @@ subroutine retrieval(files, config, fluxes, obs, states, covar)
endif
! cost state space
- ! 2Jp = (p - p0)^TB^(-1)(p - p0)
- grad_p = prod_icov(config, covar, states%px - states%px0)
- cost_p = dot_product( states%px - states%px0, real(grad_p,kind=4) )
+ if ( config%lognormal ) then
+ ! lognormal (note p is log transformed variable)
+ ! 2Jp = 2*sum(p) + (p - p0)^TB^(-1)(p - p0)
+ icovp = prod_icov(config, covar, states%px - states%px0)
+ cost_p = 2.*sum(states%px) + dot_product( states%px - states%px0, icovp )
+ else
+ ! normal
+ ! 2Jp = (p - p0)^TB^(-1)(p - p0)
+ grad_p = prod_icov(config, covar, states%px - states%px0)
+ cost_p = dot_product( states%px - states%px0, real(grad_p,kind=4) )
+ endif
write(logid,*) 'Final cost 2Jp = ',cost_p
write(logid,*) 'Final cost 2Jo = ',cost_o
diff --git a/source/simulate.f90 b/source/simulate.f90
index 71f24e2970e26b5d703ee16696b08516415571c8..d67c1a1e3aa183babc60abcafdce6a4eea3dbac5 100644
--- a/source/simulate.f90
+++ b/source/simulate.f90
@@ -408,11 +408,24 @@ subroutine simulate(iter, files, config, fluxes, obs, states, grad_o, cost_o)
if ( istate(n).eq.0 ) cycle
ns = istate(n) - minval(istate,mask=(istate.gt.0)) + 1
! state vector, px
- px((ns-1)*nbox+1:ns*nbox) = states%px((istate(n)-1)*nbox+1:istate(n)*nbox)
+ if ( config%lognormal ) then
+ ! lognormal
+ ! optimize scalar z = ln(x/xb) so x = exp(z)*xb
+ px((ns-1)*nbox+1:ns*nbox) = exp(states%px((istate(n)-1)*nbox+1:istate(n)*nbox)) * &
+ fluxes%flx_box(:,istate(n),1)
+ else
+ ! normal
+ px((ns-1)*nbox+1:ns*nbox) = states%px((istate(n)-1)*nbox+1:istate(n)*nbox)
+ endif
! transport operator, hx
hx((ns-1)*nbox+1:ns*nbox) = hx((ns-1)*nbox+1:ns*nbox) + hbox((n-1)*nbox+1:n*nbox)
end do
endif
+! if ( config%lognormal ) then
+! ! lognormal
+! ! variable is log transformed
+! px(:) = exp(px)
+! endif
endif
call system_clock(time2, countrate)
@@ -471,7 +484,13 @@ subroutine simulate(iter, files, config, fluxes, obs, states, grad_o, cost_o)
end do
ni = min(ntcini, ni)
if ( config%method.eq.'congrad'.or.config%method.eq.'m1qn3' ) then
- obs%cinipos(i) = dot_product(obs%cini(i,:), states%px(npvar+(ni-1)*ncini+1:npvar+ni*ncini))
+ if ( config%lognormal ) then
+ ! lognormal
+ obs%cinipos(i) = dot_product(obs%cini(i,:), exp(states%px(npvar+(ni-1)*ncini+1:npvar+ni*ncini)))
+ else
+ ! normal
+ obs%cinipos(i) = dot_product(obs%cini(i,:), states%px(npvar+(ni-1)*ncini+1:npvar+ni*ncini))
+ endif
else
hmat(i,npvar+(ni-1)*ncini+1:npvar+ni*ncini) = obs%cini(i,:)
obs%cinipos(i) = dot_product(hmat(i,npvar+(ni-1)*ncini+1:npvar+ni*ncini), states%px(npvar+(ni-1)*ncini+1:npvar+ni*ncini))
@@ -489,34 +508,60 @@ subroutine simulate(iter, files, config, fluxes, obs, states, grad_o, cost_o)
obs%delta(i) = obs%model(i) + obs%nee(i) + obs%fff(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%cinipos(i)
+ if ( config%lognormal ) then
+ ! optimize scalar of fluxes so no contribution from best guess fluxes (ghg)
+ obs%delta(i) = obs%model(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%cinipos(i)
+ endif
endif
! if background missing then observation should have no influence!
- if ( sum(obs%cini(i,:)).eq.0. ) then
- write(logid,*) 'WARNING: background zero for observation ',i
- obs%delta(i) = 0.
- endif
+! if ( sum(obs%cini(i,:)).eq.0. ) then
+! write(logid,*) 'WARNING: background zero for observation ',i
+! obs%delta(i) = 0.
+! endif
! calculate gradient in observation space
! ---------------------------------------
- ! Jo'(p) = H^TR^-1(H(p) - y)
- ! calculate as: Jo'(p) = sum( H_i^T*ydel_i*R_i )
call system_clock(time1, countrate)
if ( config%method.eq.'congrad'.or.config%method.eq.'m1qn3' ) then
do n = 1, ntstep
if ( istateuni(n).gt.0 ) then
- grad_o((istateuni(n)-1)*ndvar+1:istateuni(n)*ndvar) = grad_o((istateuni(n)-1)*ndvar+1:istateuni(n)*ndvar) + &
- hx((n-1)*ndvar+1:n*ndvar) * obs%delta(i) / obs%err(i)**2
+ if ( config%lognormal ) then
+ ! lognormal
+! ! Jo'(p) = diag(exp(p))H^TR^-1(H(exp(p)) - y)
+! ! calculate as: Jo'(p) = sum( exp(p)*H_i^T*ydel_i*R_i ) (note here px = exp(p))
+! grad_o((istateuni(n)-1)*ndvar+1:istateuni(n)*ndvar) = grad_o((istateuni(n)-1)*ndvar+1:istateuni(n)*ndvar) + &
+! px((n-1)*ndvar+1:n*ndvar) * hx((n-1)*ndvar+1:n*ndvar) * obs%delta(i) / obs%err(i)**2
+ ! Jo'(p) = diag(p0*exp(p))H^TR^-1(H(xb*exp(p)) - y)
+ ! calculate as: Jo'(p) = sum( xb*exp(p)*H_i^T*ydel_i*R_i ) (note here px = xb*exp(p))
+ grad_o((istateuni(n)-1)*ndvar+1:istateuni(n)*ndvar) = grad_o((istateuni(n)-1)*ndvar+1:istateuni(n)*ndvar) + &
+ px((n-1)*ndvar+1:n*ndvar) * hx((n-1)*ndvar+1:n*ndvar) * obs%delta(i) / obs%err(i)**2
+ else
+ ! normal
+ ! Jo'(p) = H^TR^-1(H(p) - y)
+ ! calculate as: Jo'(p) = sum( H_i^T*ydel_i*R_i )
+ grad_o((istateuni(n)-1)*ndvar+1:istateuni(n)*ndvar) = grad_o((istateuni(n)-1)*ndvar+1:istateuni(n)*ndvar) + &
+ hx((n-1)*ndvar+1:n*ndvar) * obs%delta(i) / obs%err(i)**2
+ endif
endif
end do
endif
+
! calculate contribution to gradient from initial mixing ratios
if ( config%opt_cini ) then
- grad_o(npvar+(ni-1)*ncini+1:npvar+ni*ncini) = grad_o(npvar+(ni-1)*ncini+1:npvar+ni*ncini) + &
- obs%cini(i,:) * obs%delta(i) / obs%err(i)**2
+ if ( config%lognormal ) then
+ ! lognormal
+ grad_o(npvar+(ni-1)*ncini+1:npvar+ni*ncini) = grad_o(npvar+(ni-1)*ncini+1:npvar+ni*ncini) + &
+ exp(states%px(npvar+(ni-1)*ncini+1:npvar+ni*ncini)) * obs%cini(i,:) * obs%delta(i) / obs%err(i)**2
+ else
+ ! normal
+ grad_o(npvar+(ni-1)*ncini+1:npvar+ni*ncini) = grad_o(npvar+(ni-1)*ncini+1:npvar+ni*ncini) + &
+ obs%cini(i,:) * obs%delta(i) / obs%err(i)**2
+ endif
endif
call system_clock(time2, countrate)
if (i.lt.10.and.iter.eq.0) write(logid,*) 'simulate: time to calc grad (ms) = ',(time2-time1)*1000/countrate
@@ -524,7 +569,7 @@ subroutine simulate(iter, files, config, fluxes, obs, states, grad_o, cost_o)
! cost observation space
! ----------------------
- ! Jo = (y - Hp)^TR^(-1)(y - Hp)
+ ! 2Jo = (y - Hp)^TR^(-1)(y - Hp)
cost_o = cost_o + obs%delta(i)**2/obs%err(i)**2
10 continue