# Domain Decomposition in space-time of 4D-VAR Data Assimilation problem: a case study on the ROMS software

Domain Decomposition of 4D-VAR Data Assimilation (DD-4DVAR) is made up of decomposition of the spate-time domain, solution of reduced forecast model and minimization of local 4D-VAR functionals. Relying on the existing software implementation of ROMS software, we describe main components of DD-4D VAR DA method, highlighting the topics that we will should address both on the mathematical problem underlying ROMS and the MPI-based code implementation of the ROMS-IS4DVAR formulation.

## Authors

• 1 publication
• 1 publication
• 2 publications
11/24/2021

### Sensitivity Analysis of Domain Decomposition method for 4D Variational Data Assimilation (DD-4DVAR DA)

We prove consistence, convergence and stability of the Domain Decomposit...
04/20/2021

### Parallel Physics-Informed Neural Networks via Domain Decomposition

We develop a distributed framework for the physics-informed neural netwo...
06/25/2020

### Domain Decomposition Parabolic Monge-Ampère Approach for Fast Generation of Adaptive Moving Meshes

A fast method is presented for adaptive moving mesh generation in multi-...
01/18/2019

### TuckerMPI: A Parallel C++/MPI Software Package for Large-scale Data Compression via the Tucker Tensor Decomposition

Our goal is compression of massive-scale grid-structured data, such as t...
02/21/2022

### Use of algebraic dual representations in domain decomposition methods for Darcy flow in 3D domain

In this work we use algebraic dual representations in conjunction with d...
03/05/2020

### Microservice Decomposition via Static and Dynamic Analysis of the Monolith

Migrating monolithic software systems into microservices requires the ap...
07/31/2019

### Balanced Identification as an Intersection of Optimization and Distributed Computing

Technology of formal quantitative estimation of the conformity of the ma...
##### This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

## 1 Roms

Regional Ocean Modeling System (ROMS) is an open-source, mature numerical framework used by both the scientific and operational communities to study ocean dynamics over 3D spatial domain and time interval.

ROMS supports different 4D-Var data assimilation (DA) methodologies such as incremental strong constraint 4D-Var (IS4D-Var) and dual formulation 4D-Var using restricted B-preconditioned Lanczos formulation of the conjugate gradient method (RBL4DVAR) [8]

. IS4D-Var and RBL4DVAR search best circulation estimate in space spanned by control vector and observations, respectively. IS4D-Var and RBL4DVAR algorithm consist of two nested loop, the outer-loop involves the module (1), namely nonlinear ROMS (NLROM) solving ROMS equations, the inner-loop involves modules (2)-(3), namely tangent linear approximation of ROMS (TLROMS) and adjoint model of ROMS (ADROMS); TLROMS and ADROMS are used for minimizing 4D-Var functional

[3] (see Figures 1-2).
NLROMS is a three-dimensional, free-surface, terrain-following ocean model that solves the Reynolds-averaged Navier-Stokes equations using the hydrostatic vertical momentum balance and Boussinesq approximation.
NLROMS computes

 xROMS(tl)=Ml−1,l(x(tl−1),f(tl),b(tl)) (1)

with the state-vector , temperature , salinity , components of vector velocity , sea surface displacement . represents nonlinear ROMS acting on , and subject to forcing , and boundary conditions during the time interval .

Minimization of the 4D-Var functional:

 JROMS(δz)=12δzB−1δz+12(Gδz−d)TR−1(Gδz−d) (2)

where are the control variable increments, is vector of innovations, , where is the observation matrix; is observation error covariance matrix and is covariance matrix of model error, is computed in the inner-loop in Figure 1.
Analysis increment, , that minimizes 4D-Var function in (2) corresponds to the solution of the equation , and is given by:

 δza=(B−1+GTR−1G)−1GTR−1d (3)

or, equivalently

 δza=BGT(GBGT+R)−1d. (4)

Equation (3) is referred to as the dual form (RBL4DVAR), while (4) is referred to as the primal form (IS4DVAR). In particular, we define

 K=BGT(GBGT+R)−1 (5)

as Kalman gain matrix.

TLROMS computes

 δxROMS(tl)≃Ml−1,lu(tl−1) (6)

where , , , and , , are the background of the circulation, surface forcing and open boundary conditions respectively, and

 u(tl−1)=((δxROMS)T(tl−1),δfT(tl),δbT(tl))T.

Equation (6) is obtained from first-order Taylor expansion of NLROMS in (1).

 u∗(tl−1)=MTl−1,lp(tl) (7)

where where is the adjoint state-vector, and are the adjoint of the surface forcing and the open boundary condition increments.

## 2 DD-4DVarDA in ROMS model

Domain Decomposition (DD) method proposed in [1] is made up of decomposition of the domain into subdomains where is the 3D spatial domain and is the time interval, solution of reduced forecast model and minimization of local 4D-Var functionals (we call it DD-4DVarDA method).
Relying on the existing software implementation, in the next we describe main components of DD-4DVarDA method, highlighting the topics that we will address both on the mathematical problem underlying ROMS and the code implementation (see steps 1-4). We focus on IS4DVAR formulation described in Section 1.

### 2.1 Decomposition in space and time of the ocean model

Decomposition of spatial domain .
We will consider a 2D decomposition of in x- and y-direction and denote the spatial domain to decompose.

ROMS uses a parallelization approach that partitions domain into tiles (see Figure 3) [2]

 Ωxy=Nsub−1⋃i=0tilei (8)

where ; and are the number of tiles set in the input file in x- and y-direction, respectively.
We denote by and the overlapping tiles regions (i.e ghost or halo area in ROMS111In ROMS, the halo area would be two grids points wide unless the MPDATA advection scheme is used, in which case it needs three.) in x- and y-direction, respectively.
Step 0: DD of in ROMS.
In our study we will assume the decomposition available in ROMS, as given in (8).

Decomposition of time interval .
ROMS does not yet implement decomposition in time direction.
Step 1: DD of in ROMS.
In our study we aim introducing a decomposition of time interval into intervals:

 Δ=Nt⋃k=1Δk:=Nt⋃k=1[t¯sk−1,t¯sk−1+Nk], (9)

where are respectively the number of subdomains of and of time such that , and .

Ocean model reduction. ROMS allows each tile (or subdomain, see Figure 3) to compute local solutions of TLROMS and ADROMS.
For and local222Let and be vectors, for simplicity of notations, we refer to as a restriction of to , i.e. and , similarly for matrix , i.e. and , according the description in [6]. TLROMS on local domain computes

 δxROMSi(tl)≃Mi,(l−1,l)ui(tl−1) (10)

and local ADROMS on local domain computes

 u∗i(tl−1)=MTi,(l−1,l)pi(tl), (11)

where , , , , are the restriction on of variables , , , and and is discrete model from to .
DD-4DVar method introduces model reduction by using the background as local initial values. For (outer loop of DD–4DVAR method [1]) do: for , posed , we let be the solution of the local model

 (12)

where , , and are respectively numbers of adjacent tiles in x- and y-direction, and are respectively background on , the vector accounting boundary conditions of and the restriction in of the matrix in (1) that is

 Mk≡M¯sr−1,¯sr:=M¯sr−1,¯sr−1+1⋯M¯sr−1,¯sr. (13)

In the following, we neglect the dependency on outer loop iteration of DD–4DVAR method. We underline that local TLROMS and ADROMS in (6) and (7) are obtained by using MPI exchange for boundary conditions, regardless of local solution on overlap area, namely they do not consider overlapping tiles conditions in (12.1) and (12.2). Consequently, we need to modify local TLROMS and ADROMS in (10) and (11) taking into account of (12.1) and (12.1) conditions. More precisely,

Step 2: TLROMS.
for , , local TLROMS on will be modified such that

 δxi,k≃Mi,kui,k−1+θI(ui,k−1)+θJ(ui,k−1) (14)

where

 θJ(ui,k−1):=∑nJiJ=1γiJ(Mi,k/HJ⋅ui,k−1/HJ−Mk/HJ⋅ui,k−1/HJ) (15)

and

 θI(ui,r−1):=∑nIiI=1γiI(Mri/HI⋅ui,r−1/HI−Mri/HI⋅ui,r−1/HI) (16)

are overlapping vectors in x- and y-direction, respectively; where parameters , denote weights.

Similarly, we need to modify local ADROMS in (11). More precisely, Step 3: ADROMS.
for , , local ADROMS in (11) on local domain will be modified such that

 ui,k−1=(Mi,r)T⋅pi,k+θJ(pi,k−1)+θI(pi,k−1) (17)

where and are defined in (15) and (16).

4D-VAR DA Operator Reduction.
Operator reduction involves the inner loop in Figure 1.
ROMS [2] computes and minimizes the operator

 JROMStilei:=JROMS/tilei (18)

where is defined in (2) and is the IS4D-Var functional in each tile with MPI exchange of boundary conditions, .

Local 4D–VAR DA functional in [1] is defined as follows

 JDD−4DVartilei×Δr:=Ji,k(xi,k)=J(xi,k)/(tilei×Δk)+OIJ(xi,k) (19)

where

 J(x)=α∥x−xb∥2B−1+∥Gx−y∥2R−1, (20)

is the DD-4DVar functional in [1],

 OIJ(xi,k)=nI∑iI=1βiI,k⋅∥xi,k/HI−xiI,k/HI∥2B−1iI+nJ∑iJ=1βiJ,k⋅∥xi,k/HJ−xiJ,k/HJ∥2B−1ij (21)

is the overlapping operator on overlapping tiles region and , and

 Ji,k(xi,k)/(tilei×Δk)=αi,k⋅∥xi,k−xMi,ki,k∥B−1i+∥Gi,rxi,k−yi,k∥2R−1i (22)

is the restriction of on , where is background, is observations vector in , is observations vector in ; , , are respectively the restrictions of covariance matrix to , and ; , are restriction of matrices and to in . Parameters , and in (21) denotes regularization parameters. We let , and .
Incremental formulation of local 4D–VAR DA functional in (19) is

 Ji,k(δzi,k)=J(δzk)/(tilei×Δk)+OIJ(δzi,k)

where are the control variable increments in .
Step 4: IS4D-Var.
From (20) local IS4D-Var function in (18) can be written

 JROMSloc=JDD−4DVartilei×Δ−OIJ (23)

Consequently, we need to add overlapping operator in order to enforce the matching of local solutions on the overlapping tiles in each time interval.

### 2.2 DD-4DVarDA in ROMS code

In Figure 6 the ROMS directory structure is shown. We focus on ROMS folder, in particular, on its folders: Tangent and Adjoint.

1. ROMS.
We need to modify routines in ROMS folder implementing decomposition as in (9) (see step 1).
Decomposition of time interval involves modification of initial conditions of nonlinear model, in TLROMS and ADROMS.
Routines involving initialization are (see Figure 6):

1. initial routine in Nonlinear folder: initializes all model variables (it is called in main3d). Note 1.1 (see step 1)
Initial routine initializes all model variable before calling main3d routine (in Nonlinear folder). Main3d routine is the main driver solving nonlinear ROMS model (background).

DD in x- and y- directions is not applied for computing background, i.e. there are not MPI communications among processes.
This means that background is computed without using DD in space, consequently we not apply DD in time.
Moreover, some values of (components velocity) are not set to zero at the end of time step, because its values on some grid points are necessary for the next time step.

2. tl_initial routine in Tangent folder: initializes all tangent model variables (it is called in i4dvar). Note 1.2 (see step 1).
We consider tl_initial routine. Parts involving initializations are:

• line 123: Initializes time stepping indices and counter.

• line 162: initialization time.

3. ad_initial routine in Adjoint folder: initializes all adjoint model variables(it is called in i4dvar). Note 1.3 (see step 1).
We consider ad_initial routine. Parts involving initializations are:

• line 113: Initializes time stepping indices and counter.

• line 152: initialization time.

Main actions to apply DD in time in ROMS are described below. Action 1.1 (see step 1).

1. Modify initials routines by adding MPI communications for initial conditions to tangent and adjoint routines in each time interval.

2. Decompose time interval (RunInterval variable in routines). One possible way forward:

• We need to use OMP threadprivate directive for replicating variables related to time interval such that each thread has its own copy. (see mod_parallel routine in Module folder at line 51 related to DD in space).
Allocate_routine in mod_parallel allocates tiles; equally we can allocate local time interval and related variables.

• We can add allocation and OMP threadprivate directive of local time intervals in mod_parallel routine. We can define first time interval (first_time_interval) and last time interval (last_time_interval) and add a for loop, after for loop involving tiles, started from first time interval up to last time interval adding a time step (dt) defined in driver.

• tl_main3d: at line 277 stars for loop involves tiles.

• ad_main3d: at line 629 stars for loop involves tiles.

• Identify tangent and adjoint variables for MPI communications in time.

3. Introduce MPI communications in time.

• We need to split ROMS MPI communicator (OCN_COMM_WORLD) for obtaining MPI communicators needed to communications among processes related to same spatial sub domain but different time intervals.
By splitting ROMS MPI communicator we obtained a new MPI communicator namely one single communicator for the current process.

Action 1.2 (see step 1).
Step 1 involves:

• tl_main3d routine: at line 142 starts while loop

on time interval (RunInterval) by increasing the step time (my_steptime). Moreover, at lines 278 and 280 it calls tl_set_massflux and tl_rho_eos, i.e. the routines we need to modify (see Action 2.1,2,3,4). Consequently, we probably need to introduce in tl_main3d MPI communications in time after the tl_set_massflux and tl_rho_eos routines.

• ad_main3d routine: at line 177 starts while loop on time interval (RunInterval) by increasing the step time (my_steptime). Moreover, at lines 632 and 634 it calls ad_rho_eos and ad_set_massflux, i.e. the routines we need to modify (see Action 3.1 and 3.2). Consequently, we probably need to introduce in ad_main3d MPI communications in time after the ad_set_massflux and ad_rho_eos routines.

2. TLROMS.
We need to modify TLROMS as in (14) (see step 2), namely we need to compute the overlapping vector in (15) and (16) and add them to tangent variables.
Action 2.1 (see step 2).
We could consider another inner loop (inside the loop over in Figure 1) over index and initial approximation of solution on local tiles at . For each iteration we need local solution on tile adjacent to each tile, this means using MPI exchange of information between adjacent tiles (in two direction see Figure 5).

We note that primitive equations of motion [5] are written in flux form transformed using orthogonal curvilinear coordinates () (see Figure 5).
We consider tl_main3d routine. In tangent folder (see Figure 6) this routine is the main driver of TLROMS configurated as a full 3D baroclinic ocean model.
tl_main3d routine calls the following subrountines.

1. tl_set_massflux calls tl_set_massflux_tile.
tl_set_massflux_tile: computes “in situ” tangent linear horizontal mass flux. Action 2.2 (see step 2).
Taking into account (15) and (16) we need to modify the code starting from line 155.

2. tl_rho_eos calls tl_rho_eos_tile.
tl_rho_eos_tile: computes “in situ” the density and other quantities (temperature, salinity,…).
Action 2.3 (see step 2).
Parts to be modified:

1. Line 505: computes ”in situ” density anomaly;

2. Line 591: computes ”in situ” Brunt-Vaisala frequency;

3. Line 1158: computes ”in situ” Brunt-Vaisala frequency;

3. tl_omega: computes vertical velocity (no modifications because DD is only horizontal, there are not overlap region along vertical direction).

Similarly to TLROMS, we need to modify ADROMS as in (17) (see step 3), namely we need to compute and add the overlapping vector in (15) and (16) to adjoint variables.
We consider ad_main3d routine. In Adjoint folder (see Figure 6) this routine is the main driver of ADROMS configurated as a full 3D baroclinic ocean model.
ad_main3d routine calls the following subrountines.

1. ad_rho_eos: computes “in situ” density and other associated quantities.
Action 3.1 (see step 3).
Parts that should be modified:

1. Lines 797 and 1758: compute ”in situ” adjoint Brut-Vaisala frequency at horizontal points.

2. Line 1070: computes ”in situ” adjoint density anomaly.

Action 3.2 (see step 3).
Part that should be modified:

1. Line 201: computes ”in situ” adjoint horizontal mass fluxes

3. ad_set_avg: accumulates and computes output time-averaged adjoint fields. (probably no modifications are needed).

4. IS4D-Var.
We need to modify routines in ROMS folder as in (23) (see step 4).
Action 4.1 (see step 4).
We need to modify IS4DVAR cost function in (2) to take in account overlap region i.e. halo region.

## 3 ROMS test cases

We have installed ROMS on the high-performance hybrid computing architecture of the Sistema Cooperativo Per Elaborazioni scientifiche multidiscipliari data center, located in the University of Naples Federico II. More precisely, the HPC architecture is made of 8 nodes, consisting of distributed memory DELL M600 blades connected by a 10 Gigabit Ethernet technology. Each blade consists of 2 Intel Xeon@2.33 GHz quadcore processors sharing 16 GB RAM memory for a total of 8 cores/blade and of 64 cores, in total.
We consider test cases referred to IS4D-Var and RBL4DVAR configured for the U.S. west coast and the California Current System (CCS). This configuration, referred to as WC13, has 30 km horizontal resolution, and 30 levels in the vertical and data assimilated every 7 days during the period July 2002 – Dec. 2004.

• IS4DVAR test cases.

To run this application I need to take the following steps.

• I customize build_roms.csh script in IS4DVAR folder and define the path to the directories where all project’s files are kept.

• I set USE_MY_LIBS to yes, then I modify source code root file Complilers/my_build_paths.sh and edit the appropriate paths for the desired compiler.

• I customize the ROMS input file (roms_wc13_2hours.in/roms_wc13_daily.in) and specify the appropriate values of tile numbers in the I-direction and J-direction i.e. NtileI=2 and NtileJ=4 (see Figure 3).

• I modify ARmake.inc in Lib folder to create the libraries for serial and parallel ARPACK.

• I customize Linux-ifort.mk (Fortran compiler is ifort) and define library locations.

• I create PBS scripts to run the test cases.

Figures 8 and 7 show results obtained applying DD in space by considering tiles and input files roms_wc13_2hours.in and roms_wc13_daily.in, respectively. They show total cost function (black curve), observation cost function (blue curve), and background cost function (red curve) and the theoretical minimum value (dashed black line) plotted on a scale. We observe that the solution has more or less converged by 25 iterations.
We will refer to the following quantities: Nouter number of outer loops, Ninner inner loop, elapsed time for each node in seconds, and maximum and minimum elapsed time among nodes in seconds.
We consider:

• case 1: Nouter=1 and Ninner=25;

• case 2: Nouter=1 and Ninner=50;

• case 3: Nouter=2 and Ninner=25;

• case 4: Nouter=2 and Ninner=50;

In Tables 1, 2 and 3, 4 we report elapsed time for each process considering roms_wc13_2hours and roms_wc13_daily.in input file, respectively.

• RBL4DVAR.

To run this application I need to take the following steps.

• I customize build_roms.csh scripts in RBL4DVAR folders and define the path to the directories where all project’s files are kept.

• I customize the ROMS input file (roms_wc13_2hours.in/roms_wc13_daily.in) and specify the appropriate values of tile numbers in the I-direction and J-direction i.e. NtileI=2 and NtileJ=4 (see Figure 3).

• I create a PBS script to run the test cases.

• IS4D-Var vs RBL4DVAR.

IS4D-Var and RBL4DVAR algorithms are based on a search for the best circulation estimate in the space spanned by the model control vector and in the dual space spanned by the observations, respectively. We note that , where is number of observations and is dimension of model control vector. Hence, the dimension of observation space is significantly smaller than the model control space, the dual formulation can reduce both memory usage and computational cost. Consequently, it appears that the dual formulation should be an easier problem to solve because of the considerably smaller dimension of the space involved, but it has practical barriers to convergence.
Consequently, we consider three applications of dual formulation:

1. using the standard preconditioning with a conjugate gradient (RBL4DVAR) method;

2. using the preconditioning and minimum residual algorithm (MINRES);

3. using a restricted B-preconditioned conjugate gradient (RPCG) approach.

Figure 9 and 10 shows the convergence of minimization algorithms RBL4DVAR and RBL4DVAR, MINRES, and RPCG compared to primal formulation IS4DVAR by considering , and , . We observe that the performance of RPCG is superior to both MINRES and RBL4DVAR, in particular, PCG ensures that RPCG converges at same rate as IS4D-Var.

• Observation impact and observation sensitivity.
The observations assimilated into the model are:

• SST satellite: Sea Surface Temperature;

• SSH satellite: Sea Surface Heights;

• Argo floats: hydrographic observations of temperature and salinity.

We consider the time average transport across 37N over the upper 500 m, denoted by , and given by

 I37N(x)=1NN∑i=1hTxi (24)

where is a vector with non-zero elements corresponding to the velocity grid points that contribute to the transport normal to the 37N section shown in Figure 12, is the number of time steps during the assimilation interval, is the model state-vector at time , and is the model time step. We consider scalar functions of ocean state vector , prior and posterior ocean state vector , where is background circulation estimate and is 4D-Var analysis. The circulation analysis increment at instant time in interval time of each 4D-Var cycle is given

 δxROMS(t)=xa(t)−xb(t), (25)

consequently, transport increment during each assimilation cycle can be expressed as:

 ΔI=Ia−Ib=I(xb+δx)−Ib≃δxT(∂J/∂x), (26)

using the tangent linear assumption in (6), (5) and (4), the circulation analysis increment is given by

 δxROMS(t)=M(t,t0)~Kd, (27)

where represents the perturbation tangent linear model for the time interval , is an approximation of defined in (5) and is innovation vector. Then, increment defined in (26) become

 ΔI≃dT~KTMT(∂I/∂x), (28)

where represents adjoint of Kalman gain matrix.
Consequently, 37N time increment of averaged transport in (24) is:

 ΔI37N≃1NdT~KTN∑l=1(Ml)Th=dTg=dT(gx+gf+gb) (29)

or

 ΔI37N≃nobs∑l=1(yl−Hl(xb(t)))gl, (30)

namely, each observation contributes to the calculation of the increase, where , contribution from initial condition increment, contribution from surface forcing increment, contribution from open boundary increment and .
Moreover, we consider the sensitivity of to variations in the observations, in particular, we consider .

Figure 10 shows the total number of observations from each observing platform that were assimilated into the model during cycle, increments and contribution of observations from each platform to

• NL: nonlinear model;

• TL: total impact;

• IC: initial condition;

• FC: surface forcing conditions;

• BC: open boundary conditions.

Figure 13 shows curve represents the evolution of the ocean state vector in time based on the either the prior or posterior control vector. In the case of the observation impact, the actual contribution of each observation, , to the change defined in (28) (blue curve vs. black curve) due to data assimilation is revealed. Conversely, the observation sensitivity quantifies the change that will occur in a (red curve vs. blue curve) as a result of perturbations in the observations, . In Figures 15, 16 and 17 we show the results obtained by running related tests on SCoPE.

## 4 Future developments

Main modification to apply DD in time in ROMS is described by point (c) in Action 1.1, i.e. introduction of MPI communications in time. More precisely, we need to introduce and manage the MPI communication in space and time. The introduction of DD in time involves the following communications among processes:

• Intra communications: by splitting MPI communicator (OCN_COMM_WORLD) to create local communicators (TASK_COMM_WORLD) to allow communications among processes related to same spatial subdomain but different time intervals.
MPI commands are

• MPI_Comm_split: partitions the group of MPI processes into disjoint subgroups and creates a new communicator (TASK_COMM_WORLD) for each subgroup.

• MPI_Isend and MPI_Irecv: sends and receives initial conditions by setting the new communicator obtained from MPI_Comm_split.

• Inter communications: by creating new communicators (OCNi_COMM_WORLD) to allow communications among processes related to different spatial subdomains but same time interval.
MPI commands are

• MPI_Intercomm_create: creates an intercommunicator for each subgroup.

• MPI_Isend and MPI_Irecv: sends and receives boundary conditions by setting the new communicator obtained from MPI_Intercomm_create.

## References

• [1] D’Amore, L., Cacciapuoti, R., 2021: “Model Reduction in Space and Time for decomposing ab initio 4D Variational Data Assimilation Problems”, Applied Numerical Mathematics, 2021, Volume 160, pp. 242-264 Elsevier, https://doi.org/10.1016/j.apnum.2020.10.003.
• [2]
• [3] Moore, A.M., H.G. Arango, G. Broquet, B.S. Powell, A.T. Weaver, and J. Zavala-Garay, 2011: The Regional Ocean Modeling System (ROMS) 4-dimensional variational data assimilation systems, Part I - System overview and formulation, Progress in Oceanography, 91, 34-49.
• [4] Hedström, K.S., 2018: Technical Manual for a Coupled Sea-Ice/Ocean Circulation Model (Version 5), OCS Study BOEM 2018-007. U.S. Department of the Interior Bureau of Ocean Energy Management Alaska OCS Region. 182 p.
• [5] Arango, H.G., Moore, A.M., Miller, A.J., Cornuelle, B.D., Lorenzo, E.Di, Cornuelle B.D. and Neilson, D.J., 2003: The ROMS Tangent Linear and Adjoint Models: A comprehensive ocean prediction and analysis system.
• [6] D’Amore, L., Arcucci, R., Pistoia, J., Toumi, R. Murli, A. 2017: On the variational data assimilation problem solving and sensitivity analysis, Journal of Computational Physics, pp. 311-326.
• [7] Langland, R. H., and N. Baker, 2004: Estimation of observation impact using the NRL atmospheric variational data assimilation adjoint system. Tellus, 56A, 189–201.
• [8] Gürol, S., Weaver, A. T., Moore, A. M., Piacentini, A., Arango, and H. G., Gratton, S. 2014: B‐preconditioned minimization algorithms for variational data assimilation with the dual formulation. Quarterly Journal of the Royal Meteorological Society, 140(679), 539-556.