Observation ensemble varies with TSMP-PDAF's domain decomposition

[jjokella]

The generated observation ensemble from PDAF generated in PDAF_enkf_obs_ensemble.F90 can vary with number of processors used in ParFlow.

The differences in the observation ensembles for different processor configurations come from possible differences in the ordering of the observations in the observation vector. F.e. for ParFlow states, the observation order given in the observation input file is re-arranged such that observation in the domain of processor 0 come first, then observation in processor 1 and so on.

On the other hand, the random numbers used to generate the ensemble always have the same order, independent of the domain decomposition or the number of processors.

Current status

For ensuring reproducible results, the user should make sure that observations are ordered in the same way, when delivered to PDAF for different domain decompositions.

Fixing this issue in the current code is not feasible due to:

1. It would need implementations in both the PDAF-library and interface.
2. Due to interface changes, the change would not be propagated into the main PDAF, but would rather stay TSMP-PDAF-specific, which should be avoided
3. We plan to move to PDAF-OMI

Plan

We resolve this issue, once PDAF-OMI is adopted in TSMP(2)-PDAF.

In PDAF-OMI, observations could receive an index as an optional attribute and this index could be used to order the random generations.

Crucial code

pdaf/src/PDAF_enkf_obs_ensemble.F90

Lines 219 to 261 in ea2b58b

	USE_OMI: IF (omi_n_obstypes == 0) THEN
	! If not using OMI
	! generate random states for local domain
	members: DO member = 1, dim_ens
	m_ens_p(:, member) = m_state_p(:)
	eigenvectors: DO j = 1, dim_obs
	CALL larnvTYPE(3, iseed, 1, randval)
	components: DO i = 1, dim_obs_p
	m_ens_p(i, member) = m_ens_p(i, member) &
	+ covar(i + local_dis(mype + 1), j) * randval
	END DO components
	END DO eigenvectors
	END DO members
	ELSE USE_OMI
	! For OMI use the matting vector map_obs_id to ensure consistency
	! if different numbers of processes are used.
	ALLOCATE(randvals(dim_obs))
	! generate random states for local domain
	membersB: DO member = 1, dim_ens
	m_ens_p(:, member) = m_state_p(:)
	! Create vector of random numbers
	DO j = 1, dim_obs
	CALL larnvTYPE(3, iseed, 1, randvals(j))
	END DO
	! Create perturbed observations
	eigenvectorsB: DO j = 1, dim_obs
	componentsB: DO i = 1, dim_obs_p
	m_ens_p(i, member) = m_ens_p(i, member) &
	+ covar(i + local_dis(mype + 1), j) * randvals(map_obs_id(j))
	END DO componentsB
	END DO eigenvectorsB
	END DO membersB
	DEALLOCATE(randvals)
	END IF USE_OMI

Main repo:
https://github.com/PDAF/PDAF/blob/642e5c0c196913bdd2a2704654441afc67421edd/src/PDAF_enkf_obs_ensemble.F90#L219-L261

Note that there is already a handling of index-ordering implmented for multiple observation types.