From cba4a79447550e71490109015d8d402eeba45156 Mon Sep 17 00:00:00 2001
From: ronesy <rlt@nilu.no>
Date: Fri, 15 Apr 2022 09:27:02 +0200
Subject: [PATCH] Update for lognormal pdf in control space and minor bug fixes
---
README.txt | 36 +----------
grid_to_ncdf/SETTINGS | 4 +-
prep_fluxes/FLUXES_ghg | 12 ++--
prep_regions/.sort.f90.swp | Bin 12288 -> 0 bytes
prep_syndata/get_error_cov.f90 | 5 +-
prep_syndata/get_flux.f90 | 33 +++++++++-
prep_syndata/get_obs.f90 | 42 +++++++-----
prep_syndata/initialize.f90 | 1 +
prep_syndata/makefile | 1 +
prep_syndata/mod_perturb.f90 | 22 +++++--
prep_syndata/mod_save.f90 | 62 ++++++++++++++----
prep_syndata/mod_settings.f90 | 15 +++++
prep_syndata/mod_var.f90 | 3 +-
prep_syndata/prep_syndata.f90 | 9 +--
settings/SETTINGS_co2_config | 11 ++--
settings/SETTINGS_co2_nest_config | 11 ++--
settings/SETTINGS_ghg_config | 17 +++--
settings/SETTINGS_ghg_nest_config | 15 +++--
source/error_cov.f90 | 25 ++++---
source/get_initcams.f90 | 104 +++++++++++++++++++++---------
source/init_cams.f90 | 41 ++++++------
source/init_ghg.f90 | 65 +++++++++++++++----
source/m1qn3_interface.f90 | 30 +++++++--
source/makefile | 8 +--
source/mod_perturb.f90 | 9 ++-
source/mod_save.f90 | 103 +++++++++++++++++++++++------
source/mod_settings.f90 | 35 +++++++---
source/read_obs.f90 | 2 +-
source/retrieval.f90 | 45 ++++++++++---
source/simulate.f90 | 75 ++++++++++++++++-----
30 files changed, 590 insertions(+), 251 deletions(-)
delete mode 100644 prep_regions/.sort.f90.swp
diff --git a/README.txt b/README.txt
index 0132ef3..fe9fe0a 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 0020b77..c4f81b0 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 f85dbb9..879ec12 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
GIT binary patch
literal 0
HcmV?d00001
literal 12288
zcmeI2Piz}S6vn5eg%&8SKv2Pj7Z9pV*6TQFC<$yUH%d~|#EERDB^;<W_Sjx$v%Adh
zCgeg72&4*56-Y=3Aq0X0To51xoIo5n@$ZDhg$o>j8-fcU@y+gfZ6{4ZNI4*8m7g}A
znfKnjZ{CbI%C3x@E6vjh{W!yQgt4RVeDlcR1N+#WON@!2?XG#g>&U*VT%K?nxw@0H
z-InEy6m}f>MCveKnO<C|=4VT_GfSuQGjr1qm&)g+D>cs!{k$i*cR}!$<vM=ecSWE#
zCq{O3sydVcN`Zf<z#%p^T8Q@2-FFSqEjMG6t4k@M6i^B%1(X6x0i}RaKq;UU_}?fH
zczf7IWV<)XxS8DdY`G^N)lDg&6i^B%1(X6x0i}RaKq;UUPzopolmbeD|DXb9gRzxs
z7`tshCXfIBcmMu><1k~F!36l_CdS?bO>hbr;EzL$eF!!|3Cx1M;G-KE`vANEP5}dq
zgHgc1pEoe}8h8~v3!VY%U=G{?_Je)k*XtSk348}W1Mh-Yz$MTEXTc=Uz)^4j{CJSD
zFTs1@C2$#B1goG53~)Qx1AaKb*caeq@G`grJg^MTf($qUe!C8H0AGVo!TaEK@H}`5
zw7@wq12W(+`2AYOz6W1{PrzH?Mern;0{4J};Ad$5E%+R~4W0u}gU5jlE`UeDL!boY
z8r=%ir4&#K?4|-kl$pxXl<Ntd3KRF=mzx-+a=6BYRX1&_@Ve!3pQ>&%*f0gB(AOE`
zPJ{e#MYv&LIh^E+ON>!;wcNm^Lh=FQtY-2%(cmW{+py9;BFhiV`f8>V-ExJjVeizL
zH9EO&1bDT5-ak?d^j^X5jzP1AB~njTGh>OIHC)hJT<_%UY|!RTe{Vcdy%p%?YOOeR
zHcOtcoZw!Kj>QcWxGwo?rfrkuL>e7S?PFsjz73D2?WPc!PS!BlIK`K2u+b2;r^6zV
z$I_39RlwK0fi`QKp&Sx9v$E+(WZKc~K)14y+wuO6*t1O1jC+r#y+oZz>Eje%b=N^L
zguS(ZUukY-rSrtRH)kx20d?8fPUdDbPrEnI&cN~*HjES7ZIvbCLbU-wh@#eER94)b
zS`($Rz-cdFl7MWB6o{M-*1^`2juIAEjvU9eL@z*p;9hGtY#J0M|2}Fr-<{(AO~Tt^
zry>rBD_d6+Vf#I%!mjr6TQxHFF#P7bWk0Zk2X;oTJ0v<Rwj2=Y`dPvlI6QixWNAlM
zI)6+b$qF8XqPL0VDb{e=Xt#=@0zpk8kJPs;(~-VUiA`#PU!pZnW$p!+WBT?0@9u$K
zPK5CAkJC&mXnjW-d~V(G7}${|ui=O)MA6|I2b3*0*L;6@K#gvSw$a@piK=)VF0v)T
zakfuRlBxAABTg^g;zG7th@*A1*i2j1Y%NhZIqMQPiG@a<Da<t<#K!Z+2`wIHjCZw5
z1ENrzEMr@~Ju;348PPtJ+ntLI4dI6tL*cnTpOh+8T@m<X(mJmPu8?Q6>DT3kc3R|$
z1{Zo{bTa;WgZb{vEXRgOe*JC=tJMx@bfj=Xql)X8RO`NyjMyj%+r`8-te{Q3mxy(u
zQj;7v(8%YU9&87JH<{0GY;5Q)C)8ch%G(x_^0}WsiF}Bb-H`mY8`=$Oo9i5vQ0LY<
zZ$OATCP8kK@J&av<t3Ws4i~7C#c%~qXfAfdqS@@+Pfduq=HgMI!2{E>@w|Yr^V5~7
zGsSYPc)B!Ksx3ocnkm)F)72`?EL5mSi^WQ<G_^EWtkB|81#e7GBeyDY!fd32l`x?5
zc!zpnsYb=}GMz70D)6`bKqLtsM2N3*Cz`vp=Ghi9!t;tSoj~R!J@iz=vj#O#(87XO
z48+pQoZGG^!<hkD0leaSdqqC)@MdUhan>31ZMvZ#*9-8<c9u+MlWbGg{W=#unxbT-
zvl)T8r%z2YZs;^jX`>Ox67~#etSfn<uLf=8bTJ3mV4(Dd2y1?K;8nCYee&dL`E8`N
zvbjC&rks|aJj_kXWT<XBsDhpzGP3D3^3gKl`3+unuksi1cKn6P^lWLNT%AmI+HOVr
E3uDnSG5`Po
diff --git a/prep_syndata/get_error_cov.f90 b/prep_syndata/get_error_cov.f90
index 4131a08..9c05c9a 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 8ec8d3a..bbd9950 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 b8661d7..3cd38b9 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 cd6cd69..42ced70 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 c1fc019..827e75a 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 e83836e..2bdd1b7 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 4b99567..84d3407 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 cf2237f..479269b 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 2aa4989..95c68f8 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 09e891b..ece7588 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 32d38c7..2b571e3 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 17d270b..f523773 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 57473da..05349ab 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 3d4b282..edbef48 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 55ac30f..d8a8886 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 aff759c..b89a396 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 0f9121f..4d6f33a 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 b4a9fd7..e6f233a 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 e814171..52a3d19 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 7c777b1..103ee50 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 133b6f6..23beee3 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 571f496..8d13924 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 baaeeed..a7942ec 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 60e054e..e1b2760 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 ffe3051..93be989 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 71f24e2..d67c1a1 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
--
GitLab