From 02f259350fd91c5c34d6db8934ea6ef3a6ff915b Mon Sep 17 00:00:00 2001
From: ronesy <rlt@nilu.no>
Date: Fri, 21 Oct 2022 13:51:23 +0200
Subject: [PATCH] updates for lognormal inversion to include influence on
mixing ratio from ocean fluxes when these are not optimized
---
source/calc_conc.f90 | 7 +++++++
source/init_ghg.f90 | 21 +++++++++++++++++++--
source/initialize.f90 | 5 +++++
source/m1qn3_interface.f90 | 6 ++++--
source/mod_fluxes.f90 | 5 +++++
source/mod_save.f90 | 13 ++++++++++---
source/simulate.f90 | 5 +++--
7 files changed, 53 insertions(+), 9 deletions(-)
diff --git a/source/calc_conc.f90 b/source/calc_conc.f90
index 807f527..e8f74af 100644
--- a/source/calc_conc.f90
+++ b/source/calc_conc.f90
@@ -101,6 +101,8 @@ subroutine calc_conc(config, fluxes, obs, ngrid, gtime, hnest, hbkg, iobs, ix1,
else
! GHG species
obs%ghg(iobs) = 0.
+ ! ocean contribution only if lognormal is true and inc_ocean is false
+ obs%ocn(iobs) = 0.
endif
do n = 1, ngrid
@@ -135,6 +137,11 @@ subroutine calc_conc(config, fluxes, obs, ngrid, gtime, hnest, hbkg, iobs, ix1,
! no nested domain
obs%ghg(iobs) = obs%ghg(iobs) + sum(hnest(:,:,n) * fluxes%flx(ix1:ix2,jy1:jy2,ilo))
endif
+ if ( config%lognormal.and..not.config%inc_ocean ) then
+ ! since for lognormal fluxes are optimized (not offsets) if ocean fluxes are not optimized need
+ ! to still account for influence on mixing ratios
+ obs%ocn(iobs) = obs%ocn(iobs) + sum(hnest(:,:,n) * fluxes%flxocn(:,:,ilo))
+ endif
endif
end do ! ngrid
diff --git a/source/init_ghg.f90 b/source/init_ghg.f90
index e6f233a..7071bde 100644
--- a/source/init_ghg.f90
+++ b/source/init_ghg.f90
@@ -83,7 +83,7 @@ subroutine init_ghg(files, config, fluxes, obs, states, covar)
integer :: yyyymmdd, hhmiss
real :: area
real(kind=8) :: jd
- integer :: ix, jy, ix0, jy0
+ integer :: ix, jy, ix0, jy0, ix1, ix2, jy1, jy2
integer :: n, nb, nread, num
integer :: ierr
real(kind=8), dimension(:), allocatable :: jdate
@@ -144,10 +144,22 @@ subroutine init_ghg(files, config, fluxes, obs, states, covar)
! flux time stamp is start of time step
fluxes%time(n) = jd
+ ! ocean fluxes for lognormal case when inc_ocean is false
+ ! cover only inversion domain
+ if ( config%lognormal.and..not.config%inc_ocean.and..not.config%nested ) then
+ ! use land-sea mask to set all land fluxes to zero
+ ix1 = int((rllx - llx + dx)/dx)
+ ix2 = int((rurx - llx)/dx)
+ jy1 = int((rlly - lly + dy)/dy)
+ jy2 = int((rury - lly)/dy)
+ fluxes%flxocn(:,:,n) = abs(lsm(:,:) - 1.)*fluxes%flx(ix1:ix2,jy1:jy2,n)
+ endif
+
jd = jd + dble(statres)
end do
+
! regional prior fluxes
! ---------------------
@@ -190,6 +202,12 @@ subroutine init_ghg(files, config, fluxes, obs, states, covar)
fluxes%flx_nest(:,:,n) = sum(flx,dim=3)/real(nread)
deallocate( flx )
+ ! ocean fluxes for lognormal case when inc_ocean is false
+ if ( config%lognormal.and..not.config%inc_ocean ) then
+ ! use land-sea mask to set all land fluxes to zero
+ fluxes%flxocn(:,:,n) = abs(lsm(:,:) - 1.)*fluxes%flx_nest(:,:,n)
+ endif
+
jd = jd + dble(statres)
end do
@@ -204,7 +222,6 @@ subroutine init_ghg(files, config, fluxes, obs, states, covar)
endif
-
! prior state vector
! ------------------
diff --git a/source/initialize.f90 b/source/initialize.f90
index 6feb513..5f682a5 100644
--- a/source/initialize.f90
+++ b/source/initialize.f90
@@ -61,6 +61,11 @@ subroutine initialize(files, config)
! -------------------
if ( config%restart.eq.1 ) then
+ inquire(file=trim(files%path_output)//trim(files%file_log),exist=lexist)
+ if ( .not.lexist ) then
+ print*, 'ERROR initialize: logfile not found for restart'
+ stop
+ endif
open(logid,file=trim(files%path_output)//trim(files%file_log),status='old',action='write',access='append',iostat=ierr)
write(logid,*) 'Picking-up run here...'
else
diff --git a/source/m1qn3_interface.f90 b/source/m1qn3_interface.f90
index 68ca430..58f6ea9 100644
--- a/source/m1qn3_interface.f90
+++ b/source/m1qn3_interface.f90
@@ -111,7 +111,7 @@ subroutine m1qn3_interface(iter, grad, files, config, fluxes, obs, states, covar
ndz=6*n+1
indic=4
dxmin=1.e-20
- df1=f
+ df1=0.01d0*f
epsrel=1.e-10
nits=maxiter
nsim=100
@@ -194,12 +194,14 @@ subroutine m1qn3_interface(iter, grad, files, config, fluxes, obs, states, covar
else
! normal
! 2Jp = x^Tx (where x = state vector in chi space)
- cost_p = real(dot_product(x,x),kind=4)
+ cost_p = real(dot_product(x,x),kind=4)
endif
! J = 0.5*(Jp + Jo)
cost = 0.5*(cost_p + cost_o)
f = cost
write(aiter,'(I2.2)') iter
+ write(logid,*) 'Cost 2Jo iteration '//aiter//' = ',cost_o
+ write(logid,*) 'Cost 2Jp iteration '//aiter//' = ',cost_p
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)
diff --git a/source/mod_fluxes.f90 b/source/mod_fluxes.f90
index fcf3918..c865b49 100644
--- a/source/mod_fluxes.f90
+++ b/source/mod_fluxes.f90
@@ -127,6 +127,11 @@ module mod_fluxes
allocate( fluxes%flx_nest(nxregrid,nyregrid,ntstate), stat = ierr )
if ( ierr.ne.0 ) stop 'ERROR alloc_fluxes: not enough memory'
endif
+ allocate( fluxes%flxocn(nxregrid,nyregrid,ntstate), stat = ierr )
+ if ( ierr.ne.0 ) then
+ write(logid,*) 'ERROR alloc_fluxes: not enough memory'
+ stop
+ endif
endif
end subroutine alloc_fluxes
diff --git a/source/mod_save.f90 b/source/mod_save.f90
index 8d13924..f905a85 100644
--- a/source/mod_save.f90
+++ b/source/mod_save.f90
@@ -92,9 +92,16 @@ module mod_save
do i = 1, nobs
cinipri = sum(obs%cini(i,:))
call caldate(obs%jdate(i), jjjjmmdd, hhmiss)
- write(100,fmt=rowfmt) obs%recs(i), jjjjmmdd, hhmiss, obs%jdate(i), &
- obs%conc(i), cinipri, obs%cinipos(i), obs%bkg(i), &
- obs%ghg(i), obs%prior(i), obs%model(i), obs%delta(i), obs%err(i)
+ if ( config%lognormal.and..not.config%inc_ocean ) then
+ ! if ocean grid cells not optimized include ocean contribution in prior and posterior model
+ write(100,fmt=rowfmt) obs%recs(i), jjjjmmdd, hhmiss, obs%jdate(i), &
+ obs%conc(i), cinipri, obs%cinipos(i), obs%bkg(i), &
+ obs%ghg(i), obs%prior(i)+obs%ocn(i), obs%model(i)+obs%ocn(i), obs%delta(i), obs%err(i)
+ else
+ write(100,fmt=rowfmt) obs%recs(i), jjjjmmdd, hhmiss, obs%jdate(i), &
+ obs%conc(i), cinipri, obs%cinipos(i), obs%bkg(i), &
+ obs%ghg(i), obs%prior(i), obs%model(i), obs%delta(i), obs%err(i)
+ endif
end do
endif
diff --git a/source/simulate.f90 b/source/simulate.f90
index d67c1a1..b31c016 100644
--- a/source/simulate.f90
+++ b/source/simulate.f90
@@ -509,8 +509,9 @@ subroutine simulate(iter, files, config, fluxes, obs, states, grad_o, cost_o)
obs%conc(i) + obs%bkg(i) + obs%cinipos(i)
else
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)
+ ! optimize scalar of fluxes so no contribution from best guess fluxes (ghg)
+ ! but add ocn contribution (zero if ocean pixels optimized)
+ obs%delta(i) = obs%model(i) - obs%conc(i) + obs%bkg(i) + obs%cinipos(i) + obs%ocn(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)
--
GitLab