Skip to content

BUG: fixed bug with multiplicative observation error #53

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
drdunmire1417 wants to merge 1 commit into
base: archive/master-kul
Choose a base branch
from
Open

BUG: fixed bug with multiplicative observation error #53

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
50 changes: 34 additions & 16 deletions lis/dataassim/algorithm/enkf/enkf_Mod.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -294,16 +294,16 @@ subroutine enkf_increments(n,k)
Nobs = Nobjs*N_obs_size
allocate(Observations(Nobs))

call generateObservations(n, k, Nobjs, Nobs, LIS_OBS_State(n,k), &
LIS_OBS_Pert_State(n,k),Observations)

!----------------------------------------------------------------------------
! Retrieve Obs_pred : model's estimate of the observations
!----------------------------------------------------------------------------

allocate(Obs_pred(Nobs,N_ens))
call LIS_surfaceModel_DAGetObsPred(n,k,Obs_pred)

call generateObservations(n, k, Nobjs, Nobs, N_ens, Obs_pred, LIS_OBS_State(n,k), &
LIS_OBS_Pert_State(n,k),Observations)

!----------------------------------------------------------------------------
! Retrieve Obs_pert : observation perturbations
!----------------------------------------------------------------------------
Expand Down Expand Up @@ -527,15 +527,13 @@ subroutine enkf_update(n,k)

Nobs = Nobjs*N_obs_size
allocate(Observations(Nobs))

call generateObservations(n, k,Nobjs, Nobs, LIS_OBS_State(n,k), &
LIS_OBS_Pert_State(n,k),Observations)


allocate(Obs_pred(Nobs,N_ens))

call LIS_surfaceModel_DAGetObsPred(n,k,Obs_pred)

call generateObservations(n, k,Nobjs, Nobs, N_ens, Obs_pred, LIS_OBS_State(n,k), &
LIS_OBS_Pert_State(n,k),Observations)


do i=1,Nobs
if(Observations(i)%assim) then
enkf_struc(n,k)%anlys_res(i) = Observations(i)%value - &
Expand Down Expand Up @@ -1371,7 +1369,7 @@ end subroutine getObsPert
! \label{generateObservations}
!
! !INTERFACE:
subroutine generateObservations(n, k, Nobjs, N_obs_size, LIS_OBS_State, &
subroutine generateObservations(n, k, Nobjs, N_obs_size, N_ens, Obs_pred, LIS_OBS_State, &
LIS_OBS_Pert_State, Observations )
! !USES:

Expand All @@ -1380,6 +1378,8 @@ subroutine generateObservations(n, k, Nobjs, N_obs_size, LIS_OBS_State, &
integer, intent(IN) :: k
integer, intent(IN) :: Nobjs
integer, intent(IN) :: N_obs_size
integer, intent(IN) :: N_ens
real, intent(in), dimension(N_obs_size,N_ens) :: Obs_pred
type(ESMF_State) :: LIS_OBS_State
type(ESMF_State) :: LIS_OBS_Pert_State
type(obs_type) :: Observations(N_obs_size)
Expand All @@ -1400,7 +1400,7 @@ subroutine generateObservations(n, k, Nobjs, N_obs_size, LIS_OBS_State, &
integer :: gval,gsize
integer :: count1
character*100 :: temp
integer :: i, t,g
integer :: i, t,g, j
character*1 :: vid(2)
integer :: status
type(ESMF_Field) :: valfield
Expand All @@ -1417,6 +1417,7 @@ subroutine generateObservations(n, k, Nobjs, N_obs_size, LIS_OBS_State, &
real :: std(N_obs_size)
integer :: species(N_obs_size)
integer :: assimflag(N_obs_size)
real :: pred_mean, pred_stddev, pred_sumsq

count1 = 1
typ = 1
Expand Down Expand Up @@ -1477,11 +1478,27 @@ subroutine generateObservations(n, k, Nobjs, N_obs_size, LIS_OBS_State, &
pertyp(count1:count1+(counts-1)) = pert_type
if(pert_type.eq.0) then
std(count1:count1+(counts-1)) = obsstd(:)
elseif(pert_type.eq.1) then

elseif(pert_type.eq.1) then ! form multiplicative observation error
do g = count1,count1+(counts-1)
if(value1(g).ne.-9999.0) then
std(g) = obsstd(g)*value1(g)
if(std(g).eq.0) std(g) = 0.1

if(std(g).eq.0) then ! if observation std = 0, calculate std of model ensembles
Copy link
Collaborator

@gdelannoy gdelannoy Oct 22, 2024 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This is a good fix for snow DA, and for our specific application. However, (1) I can imagine objections about passing information on the obs_predictions into the generation of observation specs, but do not have an efficient alternative right away. Maybe it needs to be a separate 'post'-observation routine, that kicks in for specific observation types, but I am not sure. (2) If this is done in the enkf, do we need something similar in the other DA tools? Both concerns are just to keep in mind - I am not saying that we need to do it now.

Copy link
Collaborator

@mbechtold mbechtold Oct 22, 2024

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

On (1): We could consider including in the pert_attrbs a user-defined parameter that defines a minimum obs error in the case of the multiplicative obs error option. This would make it obsolete to pass the obs_predictions. This would also solve the problem that an obs value slightly above 0 would cause a switch to the else statement which results in an abrupt change of the error from one time step to the other. This could be avoided by using this minimum obs error for low values in the observations.

Copy link
Collaborator Author

@drdunmire1417 drdunmire1417 Oct 22, 2024

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yeah, I think that Michel's suggestion is a good solution. We would probably need to pass in 2 parameters: #1 - the minimum value, below which a constant error is set, #2 - the constant error for these near-0 values

pred_mean = sum(Obs_pred(g,:)) / real(N_ens)
pred_sumsq = 0.0
do j = 1, N_ens
pred_sumsq = pred_sumsq + (Obs_pred(g,j) - pred_mean) ** 2
end do
pred_stddev = sqrt(pred_sumsq / (real(N_ens) - 1))
std(g) = pred_stddev

if(std(g).eq.0) then ! if model std also = 0, do not assimilate
std(g) = -9999.0
endif

endif

else
std(g) = -9999.0
endif
Expand All @@ -1495,6 +1512,7 @@ subroutine generateObservations(n, k, Nobjs, N_obs_size, LIS_OBS_State, &
deallocate(gid1)
deallocate(obsassimflag)
deallocate(obsstd)

enddo

!----------------------------------------------------------------------------
Expand All @@ -1515,10 +1533,10 @@ subroutine generateObservations(n, k, Nobjs, N_obs_size, LIS_OBS_State, &
Observations(t)%value = value(t)
Observations(t)%std = std(t)
Observations(t)%pert_type = pertyp(t)
if(assimflag(t).eq.1) then
Observations(t)%assim = .true.
else
if(assimflag(t).eq.0 .OR. std(t).eq.-9999.0) then
Observations(t)%assim = .false.
else
Observations(t)%assim = .true.
endif
enddo

Expand Down
10 changes: 9 additions & 1 deletion lis/dataassim/algorithm/enkf/enkf_general.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -184,7 +184,15 @@ subroutine enkf_analysis( gid, &
! compute right hand side rhs = Zpert - H*y^f for system equation,

do i=1,N_obs
rhs(i) = Observations(i)%value + Obs_err(i,n_e) - Obs_pred(i,n_e)
if(Observations(i)%pert_type.eq.0) then ! additive observation error
rhs(i) = Observations(i)%value + Obs_err(i,n_e) - Obs_pred(i,n_e)
elseif(Observations(i)%pert_type.eq.1) then ! multiplicative observation error
if(Observations(i)%value.eq.0) then ! if the observation = 0, multiply by model estimate instead
rhs(i) = Observations(i)%value + (Obs_pred_bar(i) * Obs_err(i,n_e)) - Obs_pred(i,n_e)
else
rhs(i) = Observations(i)%value + (Observations(i)%value * Obs_err(i,n_e)) - Obs_pred(i,n_e)
endif
endif
end do

! solve W*b = Zpert - H*y^f
Expand Down