The bed bottom sediment output switches in ROMS ocean.in standard input file:

! Logical switches (TRUE/FALSE) to activate writing of exposed sediment
! layer properties into HISTORY output file.  Currently, MBOTP properties
! are expected for the bottom boundary layer and/or sediment models:
!
! idBott( 1=isd50)   grain_diameter          mean grain diameter
! idBott( 2=idens)   grain_density           mean grain density
! idBott( 3=iwsed)   settling_vel            mean settling velocity
! idBott( 4=itauc)   erosion_stress          critical erosion stress
! idBott( 5=irlen)   ripple_length           ripple length
! idBott( 6=irhgt)   ripple_height           ripple height
! idBott( 7=ibwav)   bed_wave_amp            wave excursion amplitude
! idBott( 8=izdef)   Zo_def                  default bottom roughness
! idBott( 9=izapp)   Zo_app                  apparent bottom roughness
! idBott(10=izNik)   Zo_Nik                  Nikuradse bottom roughness
! idBott(11=izbio)   Zo_bio                  biological bottom roughness
! idBott(12=izbfm)   Zo_bedform              bed form bottom roughness
! idBott(13=izbld)   Zo_bedload              bed load bottom roughness
! idBott(14=izwbl)   Zo_wbl                  wave bottom roughness
! idBott(15=iactv)   active_layer_thickness  active layer thickness
! idBott(16=ishgt)   saltation               saltation height
!
!                                 1 1 1 1 1 1 1
!               1 2 3 4 5 6 7 8 9 0 1 2 3 4 5 6

Hout(idBott) == T T T T T T T T T F F F F F F F

were moved to input script sediment.in and reformatted for easy use:

! Logical switches (TRUE/FALSE) to activate writing of bed bottom sediment
! parameters, [1:Ngrids] values expected.

 Hout(isd50) == T       ! grain_diameter          mean grain diameter
 Hout(idens) == T       ! grain_density           mean grain density
 Hout(iwsed) == T       ! settling_vel            mean settling velocity
 Hout(itauc) == T       ! erosion_stress          critical erosion stress
 Hout(irlen) == T       ! ripple_length           ripple length
 Hout(irhgt) == T       ! ripple_height           ripple height
 Hout(ibwav) == T       ! bed_wave_amp            wave excursion amplitude
 Hout(izdef) == T       ! Zo_def                  default bottom roughness
 Hout(izapp) == T       ! Zo_app                  apparent bottom roughness
 Hout(izNik) == F       ! Zo_Nik                  Nikuradse bottom roughness
 Hout(izbio) == F       ! Zo_bio                  biological bottom roughness
 Hout(izbfm) == F       ! Zo_bedform              bed form bottom roughness
 Hout(izbld) == F       ! Zo_bedload              bed load bottom roughness
 Hout(izwbl) == F       ! Zo_wbl                  wave bottom roughness
 Hout(iactv) == F       ! active_layer_thickness  active layer thickness
 Hout(ishgt) == F       ! saltation               saltation height

In addition, the script sediment.in was reordered to more logical distribution of the sediment standard input parameters. However, ROMS can read the parameters in input scripts in any arbitrary order except dimension parameters at the top of ocean.in.

Similarly, the switches in input script station.in:

! Logical switches (TRUE/FALSE) to activate writing of exposed sediment
! layer properties into STATIONS output file.  Currently, MBOTP properties
! are expected for the bottom boundary layer and/or sediment models:
!
! idBott( 1=isd50)   grain_diameter          mean grain diameter
! idBott( 2=idens)   grain_density           mean grain density
! idBott( 3=iwsed)   settling_vel            mean settling velocity
! idBott( 4=itauc)   erosion_stress          critical erosion stress
! idBott( 5=irlen)   ripple_length           ripple length
! idBott( 6=irhgt)   ripple_height           ripple height
! idBott( 7=ibwav)   bed_wave_amp            wave excursion amplitude
! idBott( 8=izdef)   Zo_def                  default bottom roughness
! idBott( 9=izapp)   Zo_app                  apparent bottom roughness
! idBott(10=izNik)   Zo_Nik                  Nikuradse bottom roughness
! idBott(11=izbio)   Zo_bio                  biological bottom roughness
! idBott(12=izbfm)   Zo_bedform              bed form bottom roughness
! idBott(13=izbld)   Zo_bedload              bed load bottom roughness
! idBott(14=izwbl)   Zo_wbl                  wave bottom roughness
! idBott(15=iactv)   active_layer_thickness  active layer thickness
! idBott(16=ishgt)   saltation               saltation height
!
!                                 1 1 1 1 1 1 1
!               1 2 3 4 5 6 7 8 9 0 1 2 3 4 5 6

Sout(idBott) == T T T T T T T T T T T T T T T T

were reformatted as:

! Logical switches (TRUE/FALSE) to activate writing of exposed sediment
! layer properties into STATIONS output file.

Sout(isd50)  == T       ! grain_diameter     mean grain diameter
Sout(idens)  == T       ! grain_density      mean grain density
Sout(iwsed)  == T       ! settling_vel       mean settling velocity
Sout(itauc)  == T       ! erosion_stress     critical erosion stress
Sout(irlen)  == T       ! ripple_length      ripple length
Sout(irhgt)  == T       ! ripple_height      ripple height
Sout(ibwav)  == T       ! bed_wave_amp       wave excursion amplitude
Sout(izdef)  == T       ! Zo_def             default bottom roughness
Sout(izapp)  == T       ! Zo_app             apparent bottom roughness
Sout(izNik)  == T       ! Zo_Nik             Nikuradse bottom roughness
Sout(izbio)  == T       ! Zo_bio             biological bottom roughness
Sout(izbfm)  == T       ! Zo_bedform         bed form bottom roughness
Sout(izbld)  == T       ! Zo_bedload         bed load bottom roughness
Sout(izwbl)  == T       ! Zo_wbl             wave bottom roughness
Sout(iactv)  == T       ! active_layer_t...  active layer thickness
Sout(ishgt)  == T       ! saltation          saltation height

for easy use.

Many thanks to John Warner to suggesting this change. It makes more sense to have these output logical switches in the sediment.in input script.

WARNING: All the ocean.in, sediment.in, and station.in input scripts in the trunk and test repositories were changed. The old format for these switches will ignored in ROMS. However, if any of these variables are desired in the output NetCDF files users need the new format for these switches.

Added Staggered Grid data model (SGRID) conventions to ROMS output NetCDF files for third party post-processing tools. Check the following link for more information about ​SGRID conventions. If you read the information, you will notice that a new integer variable is defined in ROMS in spherical applications:

       int grid ;
               grid:cf_role = "grid_topology" ;
               grid:topology_dimension = 2 ;
               grid:node_dimensions = "xi_psi eta_psi" ;
               grid:face_dimensions = "xi_rho: xi_psi (padding: both) eta_rho: eta_psi (padding: both)" ;
               grid:edge1_dimensions = "xi_u: xi_psi eta_u: eta_psi (padding: both)" ;
               grid:edge2_dimensions = "xi_v: xi_psi (padding: both) eta_v: eta_psi" ;
               grid:node_coordinates = "lon_psi lat_psi" ;
               grid:face_coordinates = "lon_rho lat_rho" ;
               grid:edge1_coordinates = "lon_u lat_u" ;
               grid:edge2_coordinates = "lon_v lat_v" ;
               grid:vertical_dimensions = "s_rho: s_w (padding: none)" ; 

or in Cartesian applications:

       int grid ;
               grid:cf_role = "grid_topology" ;
               grid:topology_dimension = 2 ;
               grid:node_dimensions = "xi_psi eta_psi" ;
               grid:face_dimensions = "xi_rho: xi_psi (padding: both) eta_rho: eta_psi (padding: both)" ;
               grid:edge1_dimensions = "xi_u: xi_psi eta_u: eta_psi (padding: both)" ;
               grid:edge2_dimensions = "xi_v: xi_psi (padding: both) eta_v: eta_psi" ;
               grid:node_coordinates = "x_psi y_psi" ;
               grid:face_coordinates = "x_rho y_rho" ;
               grid:edge1_coordinates = "x_u y_u" ;
               grid:edge2_coordinates = "x_v y_v" ;
               grid:vertical_dimensions = "s_rho: s_w (padding: none)" ;

Currently, the integer value for the grid variable is arbitrary since it is not used. It is set to unity in ROMS so it can be used as a logical in some post-processing programs. Only the attributes for grid are relevant.

Two variable attributes grid and location are added to each state 2D or 3D variable. For example:

        double zeta(ocean_time, eta_rho, xi_rho) ;
                zeta:long_name = "free-surface" ;
                zeta:units = "meter" ;
                zeta:time = "ocean_time" ;
                zeta:grid = "grid" ;
                zeta:location = "face" ;
                zeta:coordinates = "lon_rho lat_rho ocean_time" ;
                zeta:field = "free-surface, scalar, series" ;
                zeta:_FillValue = 1.e+37f ;
        double u(ocean_time, s_rho, eta_u, xi_u) ;
                u:long_name = "u-momentum component" ;
                u:units = "meter second-1" ;
                u:time = "ocean_time" ;
                u:grid = "grid" ;
                u:location = "edge1" ;
                u:coordinates = "lon_u lat_u s_rho ocean_time" ;
                u:field = "u-velocity, scalar, series" ;
                u:_FillValue = 1.e+37f ;
        double v(ocean_time, s_rho, eta_v, xi_v) ;
                v:long_name = "v-momentum component" ;
                v:units = "meter second-1" ;
                v:time = "ocean_time" ;
                v:grid = "grid" ;
                v:location = "edge2" ;
                v:coordinates = "lon_v lat_v s_rho ocean_time" ;
                v:field = "v-velocity, scalar, series" ;
                v:_FillValue = 1.e+37f ;
        double temp(ocean_time, s_rho, eta_rho, xi_rho) ;
                temp:long_name = "potential temperature" ;
                temp:units = "Celsius" ;
                temp:time = "ocean_time" ;
                temp:grid = "grid" ;
                temp:location = "face" ;
                temp:coordinates = "lon_rho lat_rho s_rho ocean_time" ;
                temp:field = "temperature, scalar, series" ;
                temp:_FillValue = 1.e+37f ;

The global Conventions attribute is updated as:

               :Conventions = "CF-1.4, SGRID-0.3" ;

Many thanks to Rich Signell for bringing this to our attention and for his leading role in the standardization conventions used in post-processing tools.

This is a very important update to the ROMS 4D-Var algorithms. It includes new options and several bug corrections. We have been working and refining these algorithms for the last couple of years.

WARNING: If you update to this version, you need to recompute the error covariance normalization coefficients for your 4D-Var application. This is extremely important since several bugs were corrected in the computation of the normalization coefficients. It is the user responsibility to make sure that this takes place in order to have the correct 4D-Var solution with this ROMS update.

What is New:

• Several new C-preprocessing options are introduced:
  
  
  • BEOFS_ONLY: Use to compute the EOFs of the background error covariance.
    
    
  • BNORM: Use the background normalization of the Hessian singular vectors in various adjoint-based algorithms.
    
    
  • BGQC: Use for background quality control of observations.
    
    
  • EVOLVED_LCZ: Use to compute 4D-Var evolved Hessian singular vectors.
    
    
  • GEOPOTENTIAL_HCONV: Use to perform horizontal convolution of error covariance along geopotential surfaces in the 4D-Var algorithms.
    
    
  • LCZ_FINAL: Use to compute 4D-Var Hessian singular vectors.
    
    
  • RPCG: Use to activate the Restricted B-preconditioned Lanczos (RBLanczos) minimization method in the 4D-Var dual fourmulation algorithms (4D-PSAS and R4D-Var).
    
    
  • TIME_CONV: Use to activate weak-contraint time convolution in the 4D-Var dual formulation algorithms (4D-PSAS and R4D-Var).
    
    

• Added the option to perform a background quality control the observations in terms of the state variables index or the observation provenance for a particular state variable. It is activated with the BGQC option. Several new parameters are added to input script s4dvar.in to achieve this quality control:

  ! 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, 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: Number of observation provenances to quality control
  ! Iprovenance: Observation provenance indices to process [1:Nprovenance]
  ! P_bgqc: Standard deviation threshold [1:Nprovenance]
  
    Nprovenance == 2
    Iprovenance == 6 7 ! ARGO Temperature and Salinity
  
    P_bgqc == 4.0d0 4.0d0
  

• Add the option to include temporal decorrelation scales in the model error covariance for the 4D-Var dual formulation, weak-constraint algorithms (4D-PSAS and R4-Var). It is activated with the TIME_CONV option. Recall that both 4D-PSAS (W4DPSAS option) and R4D-Var (W4DVAR option) can be run in either strong-constraint or week-constraint modes. For strong-constraint, the user needs to set NADJ to the same value of NTIMES in ocean.in:

  ! Time-Stepping parameters.
  
    NTIMES == 192
  ...
  
  ! Output tangent linear and adjoint models parameters.
  
  ...
      NADJ == 192                              ! strong constraint
  

  Conversely in weak-constraint configurations, NADJ needs to be a multiple of NTIMES. For example:

  ! Time-Stepping parameters.
  
    NTIMES == 192
  ...
  
  ! Output tangent linear and adjoint models parameters.
  
  ...
       NADJ == 48                              ! weak constraint
  

  The parameters for the weak constraint time correlations are specified in input script s4dvar.in:

  ! 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
  

• A new 4D-Var minimization solver was added for the dual formulation (space spanned by the observation vector) drivers. Use the RPCG option to activate it in either 4D-PSAS or R4D-Var. It is named as the Restricted B-preconditioned Lanczos (RBLanczos) algorithm by Gurol et al. (2014). It is mathematically equivalent to the Restricted B-preconditioned Conjugate Gradient (RBCG) used in the primal formulation (space spanned by the model state vector). More importantly, it has similar convergence properties as the primal formulation minimization solvers making weak-constraint 4D-Var very affordable. For example, for a weak-constraint 4D-PSAS in the WC13 test application we get:
  
  
  
  The figure shows the cost functions for two outer loops with 51 inner loops each (iterations). Both outer loops are plotted continuously to facilitate comparisons. Notice that the cost functions converge rapidly around 20 interactions in the first outer loop and just 10 iterations in the second outer loop.

The RBLanczos algorithm (rpcg_lanczos.F) was very difficult to code. It gets quite complex when adding the augmented terms so can be used with multiple outer loops applications. Many thanks to Selime Gürol (CERFACS) and Andy Moore (UCSC) for helping us code this tricky solver. Please check the following paper to learn more about RBLanczos:

Gürol, S., A.T. Weaver, A.M. Moore, A. Piacentini, H.G. Arango, S. Gratton, 2014: B-preconditioned minimization algorithms for data assimilation with the dual formulation, Q. J. R. Meteorol. Soc., 140, 539-556.

• The routine normalization.F which computes the background and model error covariance (B and Q, respectively) normalization coefficient was revised and several bugs were corrected. These coefficients ensure that the diagonal elements of B and Q are equal to unity. We can use the exact (Nmethod=0; expensive) or approximated randomization (Nmethod=1; cheaper) methods. They just need to be computed once for a particular set of decorrelation scales. BEWARE, since this routine changed a lot, users need to recompute these coefficients for their 4D-Var application when updating to this version of ROMS.

• Added several new files including ad_force_dual.F, comp_Jb0.F, lanc_resid.F, rpcg_lanczos.F, sum_imp.F, time_corr.F, wrt_aug_imp.F, and wrt_evolved.F. Also, numerous updates to several subroutines were carry out.

Many thanks to Andy Moore for his help in coding and testing all these algorithms. We have been talking about this update for the last couple of years and have so many research versions of the code that were becoming unmanageable.