Skip to content
Snippets Groups Projects
Commit 4385651e authored by ronesy's avatar ronesy
Browse files

Made optimize flux offsets an option and added method M1QN3

parent 3084df20
No related branches found
No related tags found
No related merge requests found
Expand all Hide whitespace changes
Inline Side-by-side
Showing
with 1891 additions and 51 deletions
source/get_initconc.f90
+ 1
− 1
View file @ 4385651e
......@@ -117,7 +117,7 @@ subroutine get_initconc(filename, varname, concini)
! calculate geopotential height for each layer
do kz = 1, nz
geoh(kz) = (gasc*ts/gc/mmol)*log(ps/prs(kz))
geoh(kz) = (gasc*ts/gc/mmol)*log(psurf/prs(kz))
enddo
! interpolate to flexpart grid
......
source/init_ghg.f90
+ 10
− 2
View file @ 4385651e
......@@ -284,8 +284,16 @@ subroutine init_ghg(files, config, fluxes, obs, states, covar)
! allocate state and error vectors
call alloc_states(states)
! optimize offsets so prior is zero
states%px0 = 0.
! assign prior state and error vectors
if ( .not.config%offsets ) then
! optimize fluxes themselves
do n = 1, ntstate
states%px0((n-1)*nbox+1:n*nbox) = flx_box(:,n)
end do
else
! optimize flux offsets
states%px0 = 0.
endif
states%xtime(:) = fluxes%time
! average error over time steps
......
source/initialize.f90
+ 1
− 1
View file @ 4385651e
......@@ -76,7 +76,7 @@ subroutine initialize(files, config)
! check if this is a restart
! --------------------------
if ( (config%restart.eq.1).and.(config%method.eq.'congrad') ) then
if ( (config%restart.eq.1).and.(config%method.eq.'congrad'.or.config%method.eq.'m1qn3') ) then
! restart crashed run
lexist = .true.
i = -1
......
source/m1qn3_interface.f90 0 → 100644
+ 205
− 0
View file @ 4385651e
!---------------------------------------------------------------------------------------
! FLEXINVERT: m1qn3_interface
!---------------------------------------------------------------------------------------
! 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
!---------------------------------------------------------------------------------------
!
!> m1qn3_interface
!!
!! Purpose:
!! Interface with the Quasi-Newton M1QN3 algorithm
!!
!! Interface:
!! uses Reverse Protocol
!! calls simulation routine (to calculate the gradient and cost) and
!! M1QN3 iteratively until reach required precision or max no. iterations
!! n ! number of variables
!! x ! state vector, on first entry prior state vector
!! f ! value of the cost function
!! g ! gradient of the cost function
!
!---------------------------------------------------------------------------------------
subroutine m1qn3_interface(iter, grad, files, config, fluxes, obs, states, covar, cost)
use mod_var
use mod_settings
use mod_fluxes
use mod_obs
use mod_states
use mod_covar
use mod_transform
use mod_m1qn3
implicit none
type (files_t), intent (in) :: files
type (config_t), intent (in) :: config
type (fluxes_t), intent (in) :: fluxes
type (obs_t), intent (in) :: obs
type (states_t), intent (in out) :: states
type (covar_t), intent (in out) :: covar
real(kind=8), dimension(nvar,1), intent (in out) :: grad
real, intent (in out) :: cost
integer, intent (in out) :: iter
! local variables that need to be initialized
character(len=3) :: normtype
integer :: i,n,imp,io,imode(3),omode,nits,nsim,iz(5),ndz,izs(1)
integer :: indic,reverse, ierr
real :: rzs(1)
real(kind=8) :: f,dxmin,df1,epsrel,dzs(1) ! reverse protocol dzs is a dummy variable
real(kind=8), dimension(nvar) :: x, g
real(kind=8), dimension(:), allocatable :: dz
! local variables for interface with simulate.f90
real :: cost_o, cost_p
real(kind=8), dimension(nvar,1) :: phi, chi_ad, px
real, dimension(nvar) :: grad_o
character(len=2) :: aiter
! variables defined in mod_m1qn3 but need to be saved
logical :: inmemo,sscale
integer :: ntravu,id,igg,idiag,iaux,ialpha,iybar,isbar,m,mmemo,reentry,isim
real(kind=8) :: gnorm,gnorms,ps,dk,dk1
logical :: cold,warm,skip_update
integer :: itmax,moderl,jcour
real(kind=8) :: d1,t,tmin,tmax,eps1,ff,preco,precos,ys,den,ps2,hp0
! variables defined in mlis3 but need to be saved
logical :: t_increased
integer :: indica,indicd
real(kind=8) :: tesf,tesd,tg,fg,fpg,td,ta,fa,fpa,d2,fn,fp,ffn,fd
real(kind=8) :: fpd,z,test,barmin,barmul,barmax,barr,gauche,droite,taa
! initialize local variables
! --------------------------
n = nvar
! initial guess at vector that minimizes the cost in chi space
x(:) = 0d0
! initial gradient in chi space
g(:) = grad(:,1)
! initial cost
f = cost
! initialize M1QN3 settings
! -------------------------
isim=1
ndz=6*n+1
indic=4
dxmin=1.e-20
df1=f
epsrel=1.e-10
nits=maxiter
nsim=100
normtype='dfn'
io=6
imp=5
imode(1)=0
imode(2)=0
imode(3)=0
reverse=1
reentry=0
! workspace for M1QN3
allocate( dz(ndz) )
dz(:)=0d0
! interface to M1QN3
! ------------------
iterate: do
call m1qn3(n,x,f,g,dxmin,df1,&
epsrel,normtype,imp,io,imode,omode,nits,nsim,iz, &
dz,ndz,reverse,indic,isim,dzs,&
inmemo,sscale,ntravu,id,igg,idiag,iaux,ialpha,iybar,isbar,m,&
mmemo,reentry,gnorm,gnorms,ps,dk,dk1,&
cold,warm,skip_update,&
itmax,moderl,jcour,&
d1,t,tmin,tmax,eps1,ff,preco,precos,ys,den,ps2,hp0,&
t_increased,indica,indicd,&
tesf,tesd,tg,fg,fpg,td,ta,fa,fpa,d2,fn,fp,ffn,fd,&
fpd,z,test,barmin,barmul,barmax,barr,gauche,droite,taa)
iter = iter + 1
if (reverse.lt.0) exit iterate
if(iter.gt.maxiter) exit iterate
! calculate new gradient
! ----------------------
if ( (config%restart.eq.1).and.(iter.le.lastiter) ) then
! use precalculated gradient
write(aiter,fmt='(I2.2)') iter
open(100,file=trim(files%path_output)//'grad_'//aiter//'.txt',action='read',status='old',iostat=ierr)
if ( ierr.ne.0 ) then
write(logid,*) 'ERROR congrad: cannot open grad_'//aiter//'.txt'
stop
endif
do i = 1, nvar
read(100,*) grad_o(i)
end do
close(100)
! use precalculated cost obs space
open(100,file=trim(files%path_output)//'cost_obs.txt',action='read',status='old',iostat=ierr)
if ( ierr.ne.0 ) then
write(logid,*) 'ERROR congrad: cannot open cost_obs.txt'
stop
endif
do i = 1, iter
read(100,*) cost_o
end do
close(100)
! calculate cost
! Jp = x^Tx where x = state vector in chi space
cost_p = dot_product(x,x)
cost = 0.5*(cost_p + cost_o)
f = cost
else
! transform x to physical space
px = chi2phi(states, covar, x)
states%px = real(px(:,1),kind=4)
! run model
cost_o = 0.
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)
! J = 0.5*(Jp + Jo)
cost = 0.5*(cost_p + cost_o)
f = cost
write(aiter,'(I2.2)') iter
write(logid,*) 'Cost iteration '//aiter//': ',cost
endif
! transform gradient to chi space
phi(:,1) = dble(grad_o)
chi_ad = phi2chi_ad(covar, phi )
! new gradient J' = J'p + J'o
g(:) = x(:) + chi_ad(:,1)
end do iterate
end subroutine m1qn3_interface
source/makefile
+ 2
− 0
View file @ 4385651e
......@@ -42,6 +42,8 @@ SRCS = mod_var.f90 \
simulate.f90 \
congrad.f90 \
mod_analytic.f90 \
mod_m1qn3.f90 \
m1qn3_interface.f90 \
retrieval.f90 \
main.f90
......
source/mod_m1qn3.f90 0 → 100644
+ 1584
− 0
View file @ 4385651e
This diff is collapsed.
source/mod_save.f90
+ 48
− 23
View file @ 4385651e
......@@ -81,12 +81,21 @@ module mod_save
obs%nee(i), obs%fff(i), obs%ocn(i), obs%bbg(i), obs%prior(i), obs%model(i), obs%delta(i), obs%err(i)
end do
else
write(100,*) 'rec yyyymmdd hhmmss juldate conc cini bkg ghg prior post diff error'
if ( .not.config%offsets ) then
write(100,*) 'rec yyyymmdd hhmmss juldate conc cini bkg ocn prior post diff error'
else
write(100,*) 'rec yyyymmdd hhmmss juldate conc cini bkg ghg prior post diff error'
endif
do i = 1, nobs
call caldate(obs%jdate(i), jjjjmmdd, hhmiss)
write(rowfmt,'(A,I6,A)') '(A3,1X,I8,1X,I6,1X,F14.6,1X,',8,'(F10.4,1X))'
write(100,fmt=rowfmt) obs%recs(i), jjjjmmdd, hhmiss, obs%jdate(i), obs%conc(i), obs%cini(i), obs%bkg(i), &
obs%ghg(i), obs%prior(i), obs%model(i), obs%delta(i), obs%err(i)
if ( .not.config%offsets ) then
write(100,fmt=rowfmt) obs%recs(i), jjjjmmdd, hhmiss, obs%jdate(i), obs%conc(i), obs%cini(i), obs%bkg(i), &
obs%ocn(i), obs%prior(i), obs%model(i), obs%delta(i), obs%err(i)
else
write(100,fmt=rowfmt) obs%recs(i), jjjjmmdd, hhmiss, obs%jdate(i), obs%conc(i), obs%cini(i), obs%bkg(i), &
obs%ghg(i), obs%prior(i), obs%model(i), obs%delta(i), obs%err(i)
endif
end do
endif
......@@ -393,13 +402,30 @@ module mod_save
else
fpri(ix,jy,n) = fluxes%flx(ix0+ix-1,jy0+jy-1,n)/numscale
endif
! regrid and add prior flux
if ( config%nested ) then
fpos(ix,jy,n) = (states%px((n-1)*nbox+nbox_xy(ix,jy)) + fluxes%flx_nest(ix,jy,n))/numscale
if ( .not.config%offsets ) then
! optimized fluxes themselves
if ( states%px0((n-1)*nbox+nbox_xy(ix,jy)).gt.smallnum ) then
! regrid using prior flux ratio grid to box
if ( config%nested ) then
fpos(ix,jy,n) = states%px((n-1)*nbox+nbox_xy(ix,jy)) * &
fluxes%flx_nest(ix,jy,n)/states%px0((n-1)*nbox+nbox_xy(ix,jy))/numscale
else
fpos(ix,jy,n) = states%px((n-1)*nbox+nbox_xy(ix,jy)) * &
fluxes%flx(ix0+ix-1,jy0+jy-1,n)/states%px0((n-1)*nbox+nbox_xy(ix,jy))/numscale
endif
else
! don't regrid according to prior flux ratio
fpos(ix,jy,n) = states%px((n-1)*nbox+nbox_xy(ix,jy))/numscale
endif
else
fpos(ix,jy,n) = (states%px((n-1)*nbox+nbox_xy(ix,jy)) + fluxes%flx(ix0+ix-1,jy0+jy-1,n))/numscale
! optimized flux offsets so add prior fluxes
if ( config%nested ) then
fpos(ix,jy,n) = (states%px((n-1)*nbox+nbox_xy(ix,jy)) + fluxes%flx_nest(ix,jy,n))/numscale
else
fpos(ix,jy,n) = (states%px((n-1)*nbox+nbox_xy(ix,jy)) + fluxes%flx(ix0+ix-1,jy0+jy-1,n))/numscale
endif
endif
! regrid but don't add prior flux
! regrid without adjusting or 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
......@@ -411,21 +437,20 @@ module mod_save
! add ocean fluxes
! ----------------
! not needed if optimize flux offsets
! if ( nocn.gt.0 ) then
! ! ocean fluxes not optimized
! focn_reg(:,:,:) = 0.
! do n = 1, ntstate
! do jy = 1, nyregrid
! do ix = 1, nxregrid
! if ( nbox_xy(ix,jy).lt.0 ) then
! focn_reg(ix,jy,n) = fluxes%ocn_box(-1*nbox_xy(ix,jy),n)/numscale
! endif
! end do
! end do
! end do
! fpos = fpos + focn_reg
! endif
if ( .not.config%offsets.and.nocn.gt.0 ) then
! ocean fluxes not optimized
focn_reg(:,:,:) = 0.
do n = 1, ntstate
do jy = 1, nyregrid
do ix = 1, nxregrid
if ( nbox_xy(ix,jy).lt.0 ) then
focn_reg(ix,jy,n) = fluxes%ocn_box(-1*nbox_xy(ix,jy),n)/numscale
endif
end do
end do
end do
fpos = fpos + focn_reg
endif
! write ncdf file
! ---------------
......
source/mod_settings.f90
+ 5
− 5
View file @ 4385651e
......@@ -119,13 +119,13 @@ module mod_settings
integer :: datef ! end date (yyyymmdd)
character(len=3) :: spec ! species ('co2' or 'ghg')
character(len=max_name_len) :: method ! inversion method ('analytic' or 'congrad')
logical :: lognormal ! lognormal observations
logical :: offsets ! optimize offsets (true or false)
logical :: inc_ocean ! include ocean boxes in inversion (true or false)
logical :: spa_corr ! use spatial correlation (true or false)
logical :: verbose ! verbose output mode (true or false)
integer :: prior_bg ! use prior best guess (0=false, 1=true)
integer :: restart ! restart a run that crashed (0=false, 1=true)
real :: molarmass ! molar mass (in flux files, e.g. C=12, CH4=16)
real :: molarmass ! molar mass of species (corresponds to flux files, e.g. C=12, CH4=16)
real :: coeff ! conversion coefficient from ppt to ppm or ppbv
logical :: nested ! use nested flexpart output (true or false)
real :: w_edge_lon ! lon of western edge of inversion grid
......@@ -644,7 +644,7 @@ module mod_settings
logical :: match, cl
! default values for logical variables
config%lognormal = .false.
config%offsets = .false.
config%inc_ocean = .false.
config%spa_corr = .false.
config%verbose = .false.
......@@ -694,9 +694,9 @@ module mod_settings
identifier = "method:"
call read_content (line, identifier, cc, cn, cl, match)
if ( match ) config%method = cc
identifier = "lognormal:"
identifier = "offsets:"
call read_content (line, identifier, cc, cn, cl, match)
if ( match ) config%lognormal = cl
if ( match ) config%offsets = cl
identifier = "inc_ocean:"
call read_content (line, identifier, cc, cn, cl, match)
if ( match ) config%inc_ocean = cl
......
source/mod_var.f90
+ 2
− 2
View file @ 4385651e
......@@ -36,7 +36,7 @@ module mod_var
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 :: maxiter=5 ! max number of iterations for congrad
integer, parameter :: maxiter=10 ! max number of iterations for congrad
real, parameter :: trunc=1.e-15 ! truncation threshold for eigenvalues of error covariance matrix
real, parameter :: reqrd = 1.e-15 ! required reduction in gradient norm for congrad
......@@ -81,7 +81,7 @@ module mod_var
real, parameter :: rearth=6.371e6 ! radius of earth (m)
real, parameter :: mmair=28.97 ! standard molecular mass of air
real, parameter :: gasc=8.314 ! ideal gas constant (J/K/mol)
real, parameter :: ps=101325 ! surface pressure (Pa)
real, parameter :: psurf=101325 ! surface pressure (Pa)
real, parameter :: gc=9.81 ! gravitational acceleration (m/s2)
real, parameter :: ts=288.15 ! surface temperature (K)
real, parameter :: lrate=0.0065 ! lapse rate (K/m)
......
source/retrieval.f90
+ 10
− 1
View file @ 4385651e
......@@ -106,7 +106,7 @@ subroutine retrieval(files, config, fluxes, obs, states, covar)
! gradient observation space
! J'o = H^(T)R^(-1)(H(p) - y)
if ( config%method.eq.'congrad'.and.config%restart.eq.1 ) then
if ( (config%method.eq.'congrad'.or.config%method.eq.'m1qn3').and.config%restart.eq.1 ) then
open(100,file=trim(files%path_output)//'grad_00.txt',action='read',status='old',iostat=ierr)
if ( ierr.ne.0 ) then
write(logid,*) 'ERROR retrieval: cannot open grad_00.txt'
......@@ -141,6 +141,15 @@ subroutine retrieval(files, config, fluxes, obs, states, covar)
call congrad(iter, chi_ad, files, config, fluxes, obs, states, covar)
iter = 99
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
grad(:,1) = dble(grad_p + grad_o)
! transform grad to chi space
chi_ad = phi2chi_ad(covar, grad)
call m1qn3_interface(iter, chi_ad, files, config, fluxes, obs, states, covar, cost)
iter = 99
call simulate(iter, files, config, fluxes, obs, states, grad_o, cost_o)
else if ( config%method.eq.'analytic' ) then
call analytic(files, config, obs, states, covar, cost_o)
else
......
source/simulate.f90
+ 23
− 16
View file @ 4385651e
......@@ -105,7 +105,7 @@ subroutine simulate(iter, files, config, fluxes, obs, states, grad_o, cost_o)
real, dimension(:), allocatable :: px ! state variables corresponding to back trajectory time
real, dimension(:), allocatable :: hx ! transport operator at state vector resolution
real :: bkgerr, ffferr ! background and fossil fuel error in observation space
real, dimension(:,:), allocatable :: datain, dataout
real, dimension(:,:), allocatable :: datain, dataout ! used in interface to interp and average routines
! initialize
! ----------
......@@ -310,7 +310,7 @@ subroutine simulate(iter, files, config, fluxes, obs, states, grad_o, cost_o)
! CO2 tracer
! ----------
if ( config%method.eq.'congrad' ) then
if ( config%method.eq.'congrad'.or.config%method.eq.'m1qn3' ) then
! select state variables for current footprint times and
! average transport operator to state variable temporal resolution
allocate( px(ntstep*ndt*nbox), stat=ierr )
......@@ -399,7 +399,7 @@ subroutine simulate(iter, files, config, fluxes, obs, states, grad_o, cost_o)
! non-CO2 tracer
! --------------
if ( config%method.eq.'congrad' ) then
if ( config%method.eq.'congrad'.or.config%method.eq.'m1qn3' ) then
! select state variables for current footprint times and
! average transport operator to state variable temporal resolution
allocate( px(ntstep*nbox), stat=ierr )
......@@ -544,6 +544,7 @@ subroutine simulate(iter, files, config, fluxes, obs, states, grad_o, cost_o)
else
! GHG fluxes
obs%ghg(i) = 0.
obs%ocn(i) = 0.
do n = 1, ngrid
! calculate indice to fluxes
flag = 0
......@@ -560,15 +561,14 @@ subroutine simulate(iter, files, config, fluxes, obs, states, grad_o, cost_o)
obs%ghg(i) = obs%ghg(i) + sum(hnest(:,:,n) * fluxes%flx(ix1:ix2,jy1:jy2,ind))
endif
end do
! don't need separate ocean part if optimizing offsets from prior
! if ( .not.config%inc_ocean ) then
! ! account for inside domain ocean fluxes if not optimized
! obs%ocn(i) = dot_product(hocn, focn_box)
! endif
if ( .not.config%inc_ocean ) then
! account for inside domain ocean fluxes if not optimized or not lognormal
obs%ocn(i) = dot_product(hocn, focn_box)
endif
endif
! contribution from state vector
if ( config%method.eq.'congrad' ) then
if ( config%method.eq.'congrad'.or.config%method.eq.'m1qn3' ) then
obs%model(i) = dot_product(hx, px)
else
obs%model(i) = dot_product(hmat(i,:), states%px)
......@@ -590,7 +590,13 @@ 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%bbg(i) - &
obs%conc(i) + obs%bkg(i) + obs%cini(i)
else
obs%delta(i) = obs%model(i) + obs%ghg(i) - obs%conc(i) + obs%bkg(i) + obs%cini(i)
if ( .not.config%offsets ) then
! optimize fluxes themselves
obs%delta(i) = obs%model(i) + obs%ocn(i) - obs%conc(i) + obs%bkg(i) + obs%cini(i)
else
! optimize offsets so account for contribution from best guess fluxes (ghg)
obs%delta(i) = obs%model(i) + obs%ghg(i) - obs%conc(i) + obs%bkg(i) + obs%cini(i)
endif
endif
! if background missing then observation should have no influence!
......@@ -603,7 +609,7 @@ subroutine simulate(iter, files, config, fluxes, obs, states, grad_o, cost_o)
! Jo'(p) = H^TR^-1(H(p) - y)
! calculate as: Jo'(p) = sum( H_i^T*ydel_i*R_i )
if ( config%method.eq.'congrad' ) then
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) + &
......@@ -612,14 +618,10 @@ subroutine simulate(iter, files, config, fluxes, obs, states, grad_o, cost_o)
end do
endif
! if log-normal observations
! Jo'(p) = H^TW(R^(-1)(ln(y) - ln(Hp)) + 1)
! where W = 1/H
! cost observation space
! ----------------------
! Jo = (y - Hx)^TR^(-1)(y - Hx)
! Jo = (y - Hp)^TR^(-1)(y - Hp)
cost_o = cost_o + obs%delta(i)**2/obs%err(i)**2
10 continue
......@@ -663,6 +665,11 @@ subroutine simulate(iter, files, config, fluxes, obs, states, grad_o, cost_o)
write(100,fmt='(E11.4)') grad_o(i)
end do
close(100)
! save cost obs space (for m1qn3)
filename = trim(files%path_output)//'cost_obs.txt'
open(100,file=trim(filename),status='unknown',access='append',action='write',iostat=ierr)
write(100,fmt='(E11.4)') cost_o
close(100)
endif
! prepare to exit
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment