Implementing CoupledDA

[ateruzzi]

The coupled DA aims at use the physical DA (PHY DA) results to updated nutrients in biogeochemistry.
3DVar with coupled DA will:

• read the density increments provided by PHY DA
• read the pre-computed correlations between density and nutrients
• calculate the nutrient increments (nitrate and phosphate)
• update the restarts accordingly

[ateruzzi]

The name of the main subroutine has been updated from oceanvar to biovar (not necessary for coupled DA, just a matter of formality).

Starting to intorduce the reading of PHY DA increments.

[ateruzzi]

The code compiles correctly BUT it is only for testing density-nutrient DA in a simplified framework:

• The update is done on nitrate outside to the DA iteration

  3dVarBio/cnv_inn.f90

  Lines 65 to 68 in fc4b5dc

  	if (drv%densnut .eq. 1) then
  	grd%n3n(:,:,:) = 0.0
  	call ver_hor_densnut(grd%n3n)
  	endif

• Single profiles of density increments are read

  3dVarBio/get_densincr_arg.f90

  Line 1 in fc4b5dc

  subroutine get_densincr_arg

• Nitrate increments are calculated on the model grid levels

  3dVarBio/get_densincr_arg.f90

  Line 1 in fc4b5dc

  subroutine get_densincr_arg

• The horizontal operator is applied

  3dVarBio/ver_hor_densnut.f90

  Line 71 in fc4b5dc

  !********** APPLY RECURSIVE FILTERS ********** !

• Phosphate is updated accordingly and RSTafter files are written

  3dVarBio/biovar.f90

  Lines 140 to 142 in fc4b5dc

  	if (drv%dnc .eq. 1) then
  	call wrt_nut_stat
  	endif

NOTE1 that this version is for testing density-nutrent DA alone (nitrate profiles assimilation is overwritten).
NOTE2 To test the code an assimilation has to be performed (NOT nitrate)

[ateruzzi]

Compiled and tested for:

• assimilation of some O2 profiles
• apply increments on n3n and n1p based on a single density increment profile

NB: it does work only if drv%nut = 1

• density increments should be relatively dense on the vertical, as it is not interpolated by EOFs
• correlation radius quite small at the moment

[ateruzzi]

From now starting a new strategy (will overwrite the previous approach)

• The assimilation will be based on multivariate nutrient-density EOFs
• Density increments from Med-PHY will serve as pseudo-observations

The cost function will minimize the distance from nitrate observations (if availables) and from density increments (weighted by their uncertainties) and from the nitrate model background (expressed as the dot products of the EOF control space solution, i.e. the weigths v assigned to both nitrate and dens-nit EOFs.

[ateruzzi]

Makefile to be still updated and then to be compiled and tested.

[ateruzzi]

• Code has been changed to adopt the strategy of multivariate density-nutrient EOFs
• Check of consistency carried out with the help of AI https://github.com/inogs/3dVarBio/blob/CoupledDA/CODE_REVIEW_SESSION.md
• Code compiles

To be tested!

And to decide about the copyright policy

[ateruzzi]

Code tested and debugged (see https://github.com/inogs/3dVarBio/blob/CoupledDA/RECENT_CHANGES_REPORT.md):

• Tests done on assimilation density incrments only (not nitrate AND density-increment assimilation - assimilation of chl or oxy neither)
• A lot of write and print for debugging to be removed
• EOFs to be checked as well as observation error
• Copyrigth to be set

[Copilot]

Pull request overview

This PR introduces coupled data assimilation (CoupledDA) to update biogeochemical nutrients using density increments produced by a separate physical DA, by adding a new drv%dnc workflow (IO, obs-vector plumbing, EOF handling, and restart updates).

Changes:

• Add density-increment (“dnc”) observation type and vectors (dnc_str), including read/interp (get_densincr_arg / int_par_dnc) and obs operators (obs_dnc, obs_dnc_ad).
• Extend EOF/covariance handling to include density–nutrient EOFs (rdeofs_dnc, ros%neof_dnc, eva/evc_dnc) and integrate them into transforms/cost function paths.
• Rename top-level driver oceanvar → biovar and update build system inputs accordingly.

Reviewed changes

Copilot reviewed 50 out of 50 changed files in this pull request and generated 8 comments.

Show a summary per file

File	Description
wrt_nut_stat.f90	Extend nutrient write/update logic to allow drv%dnc paths.
wrt_dia.f90	Include drv%dnc conditions when writing correction outputs.
wrt_chl_stat.f90	Fix deallocation list to include LimitCorr.
ver_hor_nut_ad.f90	Route adjoint vertical EOF call to veof_dnc_ad for Var='D'.
ver_hor_nut.f90	Route vertical EOF call to veof_dnc for Var='D'.
veof_nut_ad.f90	Append dnc EOF blocks into nut EOF handling (adjoint).
veof_nut.f90	Append dnc EOF blocks into nut EOF handling (forward).
veof_dnc_ad.f90	New adjoint vertical transform for density increments.
veof_dnc.f90	New vertical transform for density increments.
tao_minimizer.f90	Adjust TAO evaluation limit and add diagnostics to drv%dia.
sav_itr.f90	Update deallocations to include grd%dnc* and ros%ev*_dnc.
resid.f90	Add density increments into obs%inc/obs%amo building.
res_inc.f90	Reset adjoint arrays for grd%dnc_ad when enabled.
readGrid.f90	Allocate grd%dnc* and read lon/lat when drv%dnc enabled.
rdeofs_multi.f90	Fix missing deallocation (x1).
rdeofs_dnc.f90	New reader for density–nutrient EOF file.
obsop_ad.f90	Call obs_dnc_ad when drv%dnc enabled.
obsop.f90	Call obs_dnc when drv%dnc enabled.
obs_vec.f90	Add dnc obs into unified observation vector.
obs_str.f90	Comment tweak.
obs_dnc_ad.f90	New adjoint obs operator for density increments.
obs_dnc.f90	New forward obs operator for density increments.
obs_arg_ad.f90	Extend N3n extension condition to include drv%dnc.
obs_arg.f90	Extend N3n extension condition to include drv%dnc.
namelists/var_3d_nml	Add neof_dnc namelist entry.
namelists/satfloat.20150101.nml	Add dnc namelist flag and set multiv=0.
mpi_str.f90	Add DncExtended_3d shared 3D extended buffer.
make3dvarwq.sh	New helper build script for 3dvarwq build.
main.f90	Switch main driver call to biovar.
int_par.f90	Call int_par_dnc when drv%dnc enabled.
grd_str.f90	Add grd%dnc and grd%dnc_ad pointers.
get_obs.f90	Call get_densincr_arg when drv%dnc enabled.
get_densincr_arg.f90	New density-increment loader + interpolation parameters (int_par_dnc).
filename_mod.f90	Add EOF_FILE_DNC and INCR_FILE defaults.
eof_str.f90	Add ros%neof_dnc and ros%ev*_dnc pointers.
drv_str.f90	Add drv%dnc flag.
dnc_str.f90	New density increment vector type/module.
def_nml_multi.f90	Read dnc from namelist and set drv%dnc.
def_nml.f90	Add neof_dnc to covariance namelist.
def_cov.f90	Read dnc EOFs and include them in ros%neof_nut when enabled.
costf.f90	Include ver_hor_nut(...,'D') path and add extra diagnostics.
cnv_inn.f90	Extend nutrient transform condition to include drv%dnc.
cnv_ctv_ad.f90	Minor formatting change.
cns_str.f90	Remove stray commented markers.
clean_mem.f90	Deallocate dnc structures and DncExtended_3d.
biovar.f90	Rename oceanvar to biovar and adjust restart-write sequencing.
RECENT_CHANGES_REPORT.md	Add generated report document.
Makefile	Add new objects and switch oceanvar.o → biovar.o.
CODE_REVIEW_SESSION.md	Add full prior review transcript.
CMakeLists.txt	Switch oceanvar.f90 → biovar.f90 in sources.

Comments suppressed due to low confidence (1)

readGrid.f90:203

• Same issue as above: (drv%dnc) is INTEGER but used as LOGICAL in .or.. Use (drv%dnc .eq. 1) to keep the condition valid Fortran.

---

💡 Add Copilot custom instructions for smarter, more guided reviews. Learn how to get started.

[Copilot]

veof_dnc is defined with no arguments, but the vertical/horizontal driver (ver_hor_nut) calls it with (NutArray, Var). As written it also writes into grd%dnc rather than the passed field, so even if the call were fixed the filtering step would still operate on an uninitialized NutArray. Please align veof_dnc with veof_nut (accept NutArray, Var and fill NutArray) or update all call sites accordingly.

[Copilot]

drv%dnc is an INTEGER flag, but here it is used as a LOGICAL term in a .or. expression. Most Fortran compilers will reject this (LOGICAL required). Please change to an explicit comparison, e.g. (drv%dnc .eq. 1).

[Copilot]

obs%no always adds dnc%nc, but when drv%dnc is disabled dnc%nc is never initialized (module variables are undefined unless explicitly set). This can yield a wrong observation count and invalid allocations. Consider obs%no = sat%nc + arg%nc + merge(dnc%nc, 0_i8, drv%dnc.eq.1) (or equivalent), or initialize dnc%nc/dnc%no to 0 unconditionally during startup.

[Copilot]

When drv%dnc==1 and drv%argo_obs==0, readGrid still allocates grd%lon/grd%lat (needed by get_densincr_arg), but clean_mem currently only deallocates grd%lon/grd%lat inside the ARGO block. Consider also deallocating grd%lon/grd%lat in the drv%dnc cleanup path (or under a combined condition) to avoid leaving them allocated.

[Copilot]

This call passes arguments to veof_dnc, but veof_dnc is currently defined with no dummy arguments and populates grd%dnc directly. This will not compile and also bypasses the NutArray that is filtered below. Recommend changing veof_dnc to take (NutArray, Var) like veof_nut, or change this call + downstream logic to use grd%dnc consistently.

[Copilot]

The OpenMP directive lists k1 as a private variable, but k1 is not declared anywhere in this subroutine. This will fail compilation with OpenMP enabled. Either declare k1 (if intended) or remove it from the PRIVATE list.

[Copilot]

These print*/write(drv%dia,*) statements run on every cost-function evaluation (potentially many times during TAO minimization). This can significantly slow the run and flood logs. Please gate them behind an existing verbosity/debug flag (e.g. drv%Verbose) or remove them before merging.

[Copilot]

Same as above: unconditional logging of ctl%f_o during each cost evaluation can create large I/O overhead. Suggest gating behind drv%Verbose/debug flag.

Suggested change

	if(MyId .eq. 0) then
	if (MyId .eq. 0 .and. drv%Verbose) then