# ROMS/TOMS Developers

Algorithms Update Web Log

arango - August 16, 2006 @ 18:02

The minimization in the descent algorithm for inner loop iteration n is given by

``` ```

where αn is the step size and dn are the conjugate direction vectors:

```  ```

or

``` ```

The steepest descent occurs when β = 0. This happens either every NiterSD iterations or when dot2 > CGtol * dot1, where

```  ```

The refined step size is computed in the second inner loop pass, m = 2, as

``` ```

arango - August 16, 2006 @ 17:59
ROMS/TOMS Incremental 4DVAR Diagram- Comments (0)

The ROMS/TOMS incremental, strong contraint 4DVAR (IS4DVAR) driver, is4dvar_ocean.h, is as follows: arango - August 16, 2006 @ 17:51
Vertical and Horizontal Convolutions- Comments (0)

Following Weaver and Courtier (2001), the vertical and horizontal correlations, C, for any state variable Φ are modeled using separated 1D and 2D diffusion equations:

```  ```

where Kv and KH are the vertical and horizontal diffusion coefficients (m2/s), assumed to be of order one. They are set to unity (isotropic correlations), if these coefficients are not found in ROMS/TOMS background/model standard deviation NetCDF file (STDname, variables Kv and Kh). The user can set Kv and Kh to the desired patterns to get anisotropic correlations. If we discretize the above diffusion operators using an explicit FTCS scheme, the stability limit is governed by:

```  ```

In practice, we set γ < 1 to be below the stability threshold, (γ = 0.5). This stability parameter is set in 4DVAR input script. There are two variables: Vgamma and Hgamma for vertical and horizontal diffusion operators. In realistic applications, the vertical difussion needs to be solved implicitly (IMPLICIT_VCONV). Otherwise, the it will be very expensive since the vertical scales are much smaller. The implicit algorithm is unconditionally stable for any time-step, so we use the maximum thickness in the above equation instead of the minimum.

The parameters controlling the vertical and horizontal length-scales of the Gaussian correlation function are:

```  ```

over a time interval 0 ≤ t ≤ T and

```  ```

where Mv and MH are the total number of vertical and horizontal integration steps. Notice that both values are forced to an even number since the squared-root diffusion operators are used to insure symmetry. The squared-root operator iterates only for half of the steps. For any given decorrelation length-scale, the number of integration steps are given by:

```  ```

The vertical and horizontal decorrelation scales (m) are set in 4DVar input script variables Vdecay(:) and Hdecay(:) for each state variable.

ROMS/TOMS reports to standard output all the parameters associated with these diffusion operators. For example,

``` 5.0000E-01  Hgamma          Horizontal diffusion stability/accuracy factor.
5.0000E-05  Vgamma          Vertical diffusion stability/accuracy factor.
20000.000  Hdecay(isFsur)  Horizontal decorrelation scale (m), free-sruface.
20000.000  Hdecay(isUbar)  Horizontal decorrelation scale (m), 2D U-momentum.
20000.000  Hdecay(isVbar)  Horizontal decorrelation scale (m), 2D V-momentum.
20000.000  Hdecay(isUvel)  Horizontal decorrelation scale (m), 3D U-momentum.
20000.000  Hdecay(isVvel)  Horizontal decorrelation scale (m), 3D V-momentum.
20000.000  Hdecay(idTvar)  Horizontal decorrelation scale (m), temp.
20000.000  Hdecay(idTvar)  Horizontal decorrelation scale (m), salt.
5.000  Vdecay(isUvel)  Vertical decorrelation scale (m), 3D U-momentum.
5.000  Vdecay(isVvel)  Vertical decorrelation scale (m), 3D V-momentum.
5.000  Vdecay(idTvar)  Vertical decorrelation scale (m), temp.
5.000  Vdecay(idTvar)  Vertical decorrelation scale (m), salt.

Minimum X-grid spacing, DXmin =  5.26835784E+00 km
Maximum X-grid spacing, DXmax =  5.58390264E+00 km
Minimum Y-grid spacing, DYmin =  5.05224079E+00 km
Maximum Y-grid spacing, DYmax =  5.35460591E+00 km
Minimum Z-grid spacing, DZmin =  1.33333333E-01 m
Maximum Z-grid spacing, DZmax =  2.86550363E+02 m

Minimum barotropic Courant Number =  2.99123135E-02
Maximum barotropic Courant Number =  8.17259140E-01
Maximum Coriolis   Courant Number =  3.43845829E-02

Horizontal convolution, NHsteps, DTsizeH =    64  3.19064213E+06 s  zeta
Horizontal convolution, NHsteps, DTsizeH =    64  3.19064213E+06 s  ubar
Horizontal convolution, NHsteps, DTsizeH =    64  3.19064213E+06 s  vbar
Horizontal convolution, NHsteps, DTsizeH =    64  3.19064213E+06 s  u
Horizontal convolution, NHsteps, DTsizeH =    64  3.19064213E+06 s  v
Horizontal convolution, NHsteps, DTsizeH =    64  3.19064213E+06 s  temp
Horizontal convolution, NHsteps, DTsizeH =    64  3.19064213E+06 s  salt

Vertical   convolution, NVsteps, DTsizeV =     6  2.05277776E+00 s  u
Vertical   convolution, NVsteps, DTsizeV =     6  2.05277776E+00 s  v
Vertical   convolution, NVsteps, DTsizeV =     6  2.05277776E+00 s  temp
Vertical   convolution, NVsteps, DTsizeV =     6  2.05277776E+00 s  salt
```

It is recommended to reduce the Vgamma parameter until a value of Mv is around six for vertical diffusion accuracy.

arango - August 16, 2006 @ 17:38
Incremental, Strong Constraint 4DVAR- Comments (0)

The incremental, strong contraint 4DVAR (IS4DVAR) approximated approach was proposed by Courtier et al. (1994) to reduce the computational cost and facilitate better preconditioning. We follow the algorithm proposed by Weaver et al (2002) that defines the increment to the difference from previous reference state.

Let xk denote the state vector on the k-th outer loop such that δxk is the increment difference from the previous state

``` ```

where xk -1 is the current estimate of the ocean state. It is equal to the background state for the first minimization (x0 = xb, δx1 = 0). The cost function on the k-th loop can be written as

```  ```

where dk – 1 = xoHxk – 1 is the innovation vector, xo is the observation vector, xb is the background state, H is a linear approximation of the observation operator H, and B and O are the background and observation error covariance matrices. Let’s introduce a new minimization variable δv, such that:

```   ```

It is more convenient to work in terms of δv (minimization space: v-space). The conditioning of the Jb in v-space is optimal since it’s Hessian, 2Jb/∂v2, yields the identity matrix. The gradient of J in v-space, denoted vJ, is given by:

``` ```

where B = BT/2B1/2 and BT/2 denotes (B1/2)T, and xJo is the gradient of Jo in model space (x-space) which is computed using the adjoint model and given by:

``` ```

Following Weaver and Courtier (2001), the background-error covariance matrix can be factored as:

``` ```

where S is a diagonal matrix of background-error standard deviations, C is a symetric matrix of background-error correlations which can be factorized as C = C1/2CT/2, G is the normalization matrix which ensures that the diagonal elements of C are equal to unity, L is a 3D self-adjoint filtering operator, and W is the grid cell area (2D fields) or volume (3D fields).

arango - August 4, 2006 @ 15:31
Background/Model Error Covariance Normalization- Comments (0)

In 4DVAR the background/model error covariance is modeled using a generalized diffusion operator. We use the approach of Weaver and Courtier (2001) to compute the correlation statistics. The normalization matrix is used to convert the covariance matrix into a correlation matrix. It insures that the diagonal elements of the background/model error covariance are equal to unity. The normalization matrix is spatially dependent and affected by the land/sea masking. Currently, there are two methods to compute these coefficients: exact (Nmethod=0) and randomization (Nmethod=1). The exact method is very expensive since the normalization coefficients are computed by perturbing each model grid cell with a delta function scaled by the area (2D) or volume (3D), and then convolving with the squared-root adjoint and tangent linear diffusion operators (ad_conv_2d.F, ad_conv_3d.F, tl_conv_2d.F, tl_conv_3d.F). The randomization (approximated) method is cheaper (Fisher and Courtier, 1995). The coefficients are initialized with random numbers having a uniform distribution (drawn from a normal distribution with zero mean and unit variance). Then, scaled by the inverse, squared-root cell area (2D) or volume (3D) and convolved with squared-root tangent linear diffusion operator. These normalization coefficients are computed in Utility/normalization.F. The correlation parameters are specified, for each state variable, in s4dvar.in. The horizontal and vertical decorrelation scales are Hdecay(:) and Vdecay(:), respectively. Check 4DVAR input script for more details. In realistic applications, it is highly recomended to do the vertical convolutions implicitly (IMPLICIT_VCONV). Otherwise, it will take too many iterations since the very different horizontal and vertical length scales.

Several plots are shown below for the normalization coefficients in the CHANNEL_NECK application. We are using a 10 km horizontal decorrelation scale and 10 m vertical decorrelation scale. This application is east-west periodic and has land/sea masking in the constriction. The normalization coefficents were computed using both exact and randomization methods. The number of iterations in the randomization method is set to Nrandom=50000 to achieve a statistical meaninfull population with approximately zero expectation mean and unit variance.

Normalization coefficients, 2D variable at RHO-points: Exact Method Normalization coefficients, 2D variable at RHO-points: Randomization Method Normalization coefficients, 2D variable at U-points: Exact Method Normalization coefficients, 2D variable at U-points: Randomization Method Normalization coefficients, 2D variable at V-points: Exact Method Normalization coefficients, 2D variable at V-points: Randomization Method Normalization coefficients, 3D variable at RHO-points, surface level: Exact Method Normalization coefficients, 3D variable at RHO-points, surface level: Randomization Method Normalization coefficients, 3D variable at U-points, surface level: Exact Method Normalization coefficients, 3D variable at U-points, surface level: Randomization Method Normalization coefficients, 3D variable at V-points, surface level: Exact Method Normalization coefficients, 3D variable at V-points, surface level: Randomization Method 