s4dvar.in

4DVar Model Input Script - s4dvar.in

The s4dvar.in file sets the parameters for the ROMS 4-Dimensional Variational Data Assimilation (4DVar) algorithms. The name of this file is set by the APARNAM keyword in the ocean.in file. This standard input ASCII file is organized in several sections as shown below, with links to more detailed explanation where required.

Notice: Detailed information about ROMS input script file syntax can be found here.

Notice: A default s4dvar.in input script is provided in the User/External subdirectory.

General Controls and Parameters

• Number of iterations in the biconjugate gradient algorithm used to solve the elliptic equation for sea surface height in the error covariance balance operator: [1:Ngrids] values expected.

Nbico == 200

• Parameters used to compute balanced salinity in terms of temperature using empirical T-S relationships in the error covariance balance operator: [1:Ngrids] values expected.

dTdz_min == 0.001d0  ! minimun dT/dz (Celsius/m)
ml_depth == 100.0d0  ! mixed-layer depth (m; positive)

• Balance operator level of no motion depth (m) used when computing balanced free-surface contribution: [1:Ngrids] values expected.

LNM_depth == 1000.0d0  ! meters, positive

• Balance operator level of no motion flag used to compute balanced free-surface contribution:

! [0] Integrate from local bottom to the surface
! [1] Integrate from LNM_depth to surface or integrate from local bottom
! if shallower than LNM_depth

LNM_flag = 0

• Balance operator logical switches for state variables to consider in the error covariance multivariate constraints.

balance(isSalt) = T  ! salinity
balance(isFsur) = T  ! free-sruface
balance(isVbar) = F  ! 2D momentum (ubar, vbar)
balance(isVvel) = T  ! 3D momentum (u, v)

• Parameter to process the Nvct eigenvector of the stabilized representer matrix when computing array modes (here, Nvct=Ninner is the most important while Nvct=1 is the least important) OR cut-off parameter for the clipped analysis to disregard potentially unphysical array modes (that is, all the eigenvectors < Nvct are disregarded).

Nvct = 50

• Upper bound on the relative error of the gradient for the Lanczos conjugate gradient algorithm.

GradErr = 1.0d-4

• Maximum error bound on Hessian eigenvectors in the Lanczos conjugate gradient algorithm. Note that even quite inaccurate eigenvectors are useful for pre-conditioning purposes.

HevecErr = 1.0d-1

• Switch (T/F) to compute approximated Hessian eigenpairs in the Lanzos conjugate gradient algorithm.

LhessianEV = T

• Switch (T/F) to activate hot start in weak-constraint (W4DVAR and W4DPSAS) algorithms of subsequent outer loops.

LhotStart = T

• Switch (T/F) to activate IS4DVAR conjugate gradient preconditioning. Two types of Limited-Memory Preconditioner (LMP) are available (Tshimanga et al., 2008): spectral LMP and Ritz LMP.

Lprecond = F

• Switch to activate either Ritz Limited-Memory Preconditioner (T) or spectral Limited-Memory Preconditioner (F) to the IS4DVAR algorithm.

Lritz = T

• If preconditioning, specify number of eigenpairs to use. If zero, use HevecErr parameter to determine the number of converged eigenpairs.

NritzEV = 0

• If weak constraint 4DVar, set number of iterations in the Lanczos algorithm used to estimate the posterior analysis error covariance matrix.

NpostI = 50

• If observations impact or observations sensitivity, set the 4D-Var outer loop to process. It must be less than or equal to Nouter.

Nimpact = 1

• Number of extra-observation classes (NextraObs), observation type indices (ExtraIndex), and observation type names (ExtraName) to consider in addition to the 1-to-1 associated with the state variables. It is used in observation operators that require more that one state variable to evaluate a particular observation type like HF radials, travel time, pressure, etc.
  
  In any application, the number of observation types is computed as: NobsVar = NstateVar + NextraObs.
  
  NextraObs values are expected for keywords ExtraIndex and ExtraName. If NextraObs > 1, enter one observation type name per line and use a backslash as the continuation.

NextraObs = 0
ExtraIndex = 20
ExtraName = radial

• If weak constraint 4DVar, set diffusive relaxation coefficients (m²/s) used to relax representer tangent linear solution to privious Picard iteration linearized trajectory.

tl_M2diff == 0.0d0  ! 2D momentum
tl_M3diff == 100.0d0  ! 3D momentum

tl_Tdiff == 50.0d0 50.0d0  ! NT tracers

Normalization Options and Parameters

• Switches (T/F) to create and write error covariance normalization factors for model, initial conditions, boundary conditions, and surface forcing. If TRUE, these factors are computed and written to NRMname(1:4) NetCDF files. If FALSE, they are read from NRMname(1:4) NetCDF file. The computation of these factors is very expensive and need to be computed only once for a particular application provided that grid land/sea masking, and decorrelation scales remains the same. Notice that four values are needed (1=initial conditions, 2=model, 3=boundary conditions, 4=surface forcing) per each nested grid: [1:4,1:Ngrids] values expected.

LdefNRM == F F F F  ! Create a new normalization files
LwrtNRM == F F F F  ! Compute and write normalization

• Switches to compute the correlation normalization coefficients for model error covariance.

CnormM(isFsur) = T  ! 2D variable at RHO-points
CnormM(isUbar) = T  ! 2D variable at U-points
CnormM(isVbar) = T  ! 2D variable at V-points
CnormM(isUvel) = T  ! 3D variable at U-points
CnormM(isVvel) = T  ! 3D variable at V-points
CnormM(isTvar) = T T  ! NT tracers

• Switches to compute the correlation normalization coefficients for initial conditions error covariance.

CnormI(isFsur) = T  ! 2D variable at RHO-points
CnormI(isUbar) = T  ! 2D variable at U-points
CnormI(isVbar) = T  ! 2D variable at V-points
CnormI(isUvel) = T  ! 3D variable at U-points
CnormI(isVvel) = T  ! 3D variable at V-points
CnormI(isTvar) = T T  ! NT tracers

• Switches to compute the correlation normalization coefficients for boundary conditions error covariance.

CnormB(isFsur) = T  ! 2D variable at RHO-points
CnormB(isUbar) = T  ! 2D variable at U-points
CnormB(isVbar) = T  ! 2D variable at V-points
CnormB(isUvel) = T  ! 3D variable at U-points
CnormB(isVvel) = T  ! 3D variable at V-points
CnormB(isTvar) = T T  ! NT tracers

• Switches to compute the correlation normalization coefficients for surface forcing error covariance.

CnormF(isUstr) = T  ! surface U-momentum stress
CnormF(isVstr) = T  ! surface V-momentum stress
CnormF(isTsur) = T T  ! NT surface tracers flux

• Correlation normalization method:

! [0] Exact, very expensive
! [1] Approximated, randomization

Nmethod == 0

• If randomization, select random number generation scheme:

! [1] Gaussian distributed deviates, numerical recipes

Rscheme == 1

• Number of iterations to compute correlation normalization coefficients via the randomization approach. A large number is required to be statistically meaningful and achieve zero expectation mean and unit variance. These factors insure that the error covariance diagonal elements are equal to unity.

Nrandom = 5000

Error Covariance Parameters

• Horizontal and vertical stability and accuracy factors (< 1) used to time-step discretized convolution operators below their theoretical limit. Notice that four values [1:4] are needed for each factor to facilitate the error covariance modeling: 1=initial conditions, 2=model, 3=boundary conditions, and 4=surface forcing.

! IC Model OBC Sur For

Hgamma = 0.5 0.5 0.5 0.5  ! horizontal operator
Vgamma = 0.05 0.05 0.05 0.05  ! vertical operator

• Model error covariance: horizontal, isotropic decorrelation scales (m). These scales are only used in weak-constraint data assimilation.

HdecayM(isFsur) == 50.0d+3  ! free-surface
HdecayM(isUbar) == 50.0d+3  ! 2D U-momentum
HdecayM(isVbar) == 50.0d+3  ! 2D V-momentum
HdecayM(isUvel) == 50.0d+3  ! 3D U-momentum
HdecayM(isVvel) == 50.0d+3  ! 3D V-momentum
HdecayM(isTvar) == 50.0d+3 50.0d+3  ! 1:NT tracers

• Model error covariance: vertical, isotropic decorrelation scales (m).

VdecayM(isUvel) == 100.0d0  ! 3D U-momentum
VdecayM(isVvel) == 100.0d0  ! 3D V-momentum
VdecayM(isTvar) == 100.0d0 100.0d0  ! 1:NT tracers

• Model error covariance: temporal decorrelation scales (days). This scales are only used in weak-constraint data assimilation.

TdecayM(isFsur) == 1.0d0  ! free-surface
TdecayM(isUbar) == 1.0d0  ! 2D U-momentum
TdecayM(isVbar) == 1.0d0  ! 2D V-momentum
TdecayM(isUvel) == 1.0d0  ! 3D U-momentum
TdecayM(isVvel) == 1.0d0  ! 3D V-momentum
TdecayM(isTvar) == 1.0d0 1.0d0  ! 1:NT tracers

• Initial conditions error covariance: horizontal, isotropic decorrelation scales (m).

HdecayI(isFsur) == 100.0d+3  ! free-surface
HdecayI(isUbar) == 100.0d+3  ! 2D U-momentum
HdecayI(isVbar) == 100.0d+3  ! 2D V-momentum
HdecayI(isUvel) == 100.0d+3  ! 3D U-momentum
HdecayI(isVvel) == 100.0d+3  ! 3D V-momentum
HdecayI(isTvar) == 100.0d+3 100.0d+3  ! 1:NT tracers

• Initial conditions error covariance: vertical, isotropic decorrelation scales (m).

VdecayI(isUvel) == 100.0d0  ! 3D U-momentum
VdecayI(isVvel) == 100.0d0  ! 3D V-momentum
VdecayI(isTvar) == 100.0d0 100.0d0  ! 1:NT tracers

• Boundary conditions error covariance: horizontal, isotropic decorrelation scales (m). A value is expected for each boundary edge in the following order:

! 1: west 2: south 3: east 4: north

HdecayB(isFsur) == 100.0d+3 100.0d+3 100.0d+3 100.0d+3  ! free-surface
HdecayB(isUbar) == 100.0d+3 100.0d+3 100.0d+3 100.0d+3  ! 2D U-momentum
HdecayB(isVbar) == 100.0d+3 100.0d+3 100.0d+3 100.0d+3  ! 2D V-momentum
HdecayB(isUvel) == 100.0d+3 100.0d+3 100.0d+3 100.0d+3  ! 3D U-momentum
HdecayB(isVvel) == 100.0d+3 100.0d+3 100.0d+3 100.0d+3  ! 3D V-momentum
HdecayB(isTvar) == 4*100.0d+3 4*100.0d+3  ! 1:NT tracers

• Boundary conditions error covariance: vertical, isotropic decorrelation scales (m). A value is expected for each boundary edge in the following order:

! 1: west 2: south 3: east 4: north

VdecayB(isUvel) == 100.0d0 100.0d0 100.0d0 100.0d0  ! 3D U-momentum
VdecayB(isVvel) == 100.0d0 100.0d0 100.0d0 100.0d0  ! 3D V-momentum
VdecayB(isTvar) == 4*100.d0 4*100.d0  ! 1:NT tracers

• Surface forcing error covariance: horizontal, isotropic decorrelation scales (m).

HdecayF(isUstr) == 100.0d+3  ! surface U-momentum stres
HdecayF(isVstr) == 100.0d+3  ! surface V-momentum stress
HdecayF(isTsur) == 100.0d+3 100.0d+3  ! 1:NT surface tracers flux

• Select flag for BackGround Quality Control (BGQC) of observations:
  
  
  [1] Quality control in terms of state variable indices
  [2] Quality control in terms of observation provenance

bgqc_type == 1

• If BGQC is in terms of state variables, set number of standard deviations threshold to use in the quality control rejection of observations.
  
  
  Use a large value (say, 1.0d+5) if you do not want to reject observations associated with a particular state variable.
  
  Use a small value (typically, 4 standard deviations) to perform quality control of observations for a particular state variable.

S_bgqc(isFsur) == 1.0d+5  ! free-surface
S_bgqc(isUbar) == 1.0d+5  ! 2D U-momentum
S_bgqc(isVbar) == 1.0d+5  ! 2D V-momentum
S_bgqc(isUvel) == 1.0d+5  ! 3D U-momentum
S_bgqc(isVvel) == 1.0d+5  ! 3D V-momentum
S_bgqc(isTvar) == 4.0d0 4.0d0  ! 1:NT tracers

• If BGQC is in terms of observation provenance (OP), set number of standard deviations threshold to use in the quality control rejection of observations.
  
  
  Use a small value (say, 4 standard deviations) to perform quality control for the desired observation provenance(s).

Nprovenance == 2  ! # of OP to quality control
Iprovenance == 6 7  ! OP indicies to process [1:Nprovenance]
P_bgqc == 4.0d0 4.0d0  ! Standard deviation threshold [1:Nprovenance]

• If applicable, set switches (T/F) used to adjust surface tracer flux: [1:NT,1:Ngrids] values expected.

Lstflux == T T  ! NT tracers

Boundary Condition Parameters

If applicable, set switches to adjust state variables at the open boundaries. Notice that a value is expected for each boundary segment per nested grid: [1:4,1:Ngrids] values expected. The boundary order is: 1=west, 2=south, 3=east, and 4=north. That is, anticlockwise starting at the western boundary.

• When processing momentum, you need to activate both components. If processing 2D momentum, you need to activate both free-surface and 3D-momentum at the processing boundary.

! W S E N ____N____
! e o a o | 4 |
! s u s r | |
! t t t t 1 W E 3
! h h | |
! |____S____|
! 1 2 3 4 2

Lobc(isFsur) == F F F F  ! free-surface
Lobc(isUbar) == F F F F  ! 2D U-momentum
Lobc(isVbar) == F F F F  ! 2D V-momentum
Lobc(isUvel) == F F F F  ! 3D U-momentum
Lobc(isVvel) == F F F F  ! 3D V-momentum

• If applicable, set switches to adjust state tracer variables at the open boundaries. Notice that a value is expected for each tracer at each boundary segment per nested grid: [1:4,1:NT,1:Ngrids] values expected. The boundary order is the same as above. Notice that the first line has the values for temperature boundaries, the second is salinity, and so on.

Lobc(isTvar) == F F F F \
F F F F

Input and Output NetCDF Files

• Input model, initial conditions, boundary conditions, and surface forcing standard deviation file names: [1:Ngrids] values expected.

STDnameM == ocean_std_m.nc
STDnameI == ocean_std_i.nc
STDnameB == ocean_std_b.nc
STDnameF == ocean_std_f.nc

• Input/output model, initial conditions, boundary conditions, and surface forcing error covariance normalization factors file name: [1:Ngrids] values expected.

NRMnameM == ocean_nrm_m.nc
NRMnameI == ocean_nrm_i.nc
NRMnameB == ocean_nrm_b.nc
NRMnameF == ocean_nrm_f.nc

• Input/output observation file name: [1:Ngrids] values expected.

OBSname == ocean_obs.nc

• Input/output Hessian eigenvectors file name: [1:Ngrids] values expected.

HSSname == ocean_hss.nc

• Input/output Lanczos vectors file name: [1:Ngrids] values expected.

LCZname == ocean_lcz.nc

• Output time-evolved Lanczos vectors file name: [1:Ngrids] values expected.

LZEname == ocean_lze.nc

• Output model data at observation locations file name: [1:Ngrids] values expected.

MODname == ocean_mod.nc

• Output posterior error covariance matrix file name: [1:Ngrids] values expected.

ERRname == ocean_err.nc