ROMS/TOMS Developers

Added a new incremental, strong contraint 4DVAR driver which uses the conjugate gradient algorithm proposed by Mike Fisher (ECMWF). Many thanks to Mike for providing us his unpublished manuscript (Efficient Minimization of Quadratic Penalty Functions, 1997) and prototype source code. Many thanks to Andy Moore for his help in the debugging of the new algorithm.

The new driver is more efficient since only one call to the tangent linear model is needed per inner loop iteration. The old algorithm required a second call to compute the optimum conjugate gradient algorithm step size. The new algorithm can be activated with cpp option IS4DVAR and the old with IS4DVAR_OLD.

The developing of this new driver has been done in two stages:

• Basic conjugate gradient algorithm where the new initial conditions increment and its associated gradient are adjusted to the optimum line search. In addition, the new gradient at each inner loop iteration is orthogonalized against all previous gradient using the modified Gramm-Schmidt algorithm.
• Lanczos algorithm to compute the eigenvectors and eigenvalues used to estimate the minimization Hessian matrix for preconditioning.

Currently, only the first stage is available. The new conjugate gradient cgradient.F follows Fisher (1997) algorithm:

Given an initial increment X(0), gradient G(0), descent direction d(0)=-G(0), and trial step size τ(0) , the minimization algorithm at the k-iteration is:

• Run tangent linear model initialized with trial step, Xhat(k), Eq 5a:

  Xhat(k) = X(k) + tau(k) * d(k)

• Run adjoint model to compute gradient at trial point, Ghat(k), Eq 5b:

  Ghat(k) = GRAD[ f(Xhat(k)) ]

• Compute optimum step size, α(k), Eq 5c:

  α(k) = τ(k) * <d(k),G(k)> / (<d(k),G(k)> – <d(k),Ghat(k)>)

  here <…> denotes dot product.

• Compute new starting point (TLM increments), X(k+1), Eq 5d:

  X(k+1) = X(k) + α(k) * d(k)

• Compute gradient at new point, G(k+1), Eq 5e:

  G(k+1) = G(k) + (alpha(k) / tau(k)) * (Ghat(k) – G(k))

• Orthogonalize new gradient, G(k+1), against all previous gradients [G(k), …, G(0)], in reverse order, using the modified Gramm-Schmidt algorithm. Notice that we need to save all inner loop gradient solutions. Currently, all the gradients are stored in the adjoint NetCDF file(s) and not in memory.
• Compute new descent direction, d(k+1), Eqs 5g and 5f:

  β(k+1) = <G(k+1),G(k+1)> / <G(k),G(k)>

  d(k+1) = – G(k+1) + β(k+1) * d(k)

  After the first iteration, the trial step size is:

  τ(k) = α(k-1)

The new driver is is4dvar_ocean.h and the old becomes is4dvar_ocean_old.h. The old driver is kept for comparison and will become obsolete in the future. A new generic dot product routine state_dotprod.F is introduced.

The new algorithm is still under testing but the results are promising as shown below for the double gyre test case (outer=1, inner=80 for both new and old drivers). Notice that the gradient norm in the new algorithm decreased by about five-orders of magnitude compared to only two-orders of magnitude in the old algorithm. The new algorithm is much cheaper.

The minimization of the cost function is also improved.

Added a VERIFICATION option that allows to compute nonlinear model solution at observation locations. This option is independent of any of the variational data assimilation options. However, it uses the same observation interface. Like in the 4DVAR drivers, the model solution at the observation locations and the model-observation statistics are written to MODname NetCDF file. Many thanks to Brian Powell for this suggestion.

I modified and updated several routines:

• Modified nl_ocean.h to allow computing and saving of the nonlinear model solution at the observation spatial and temporal locations. Several routines associated with data assimilation were changed slightly to accomplish this: cppdefs.h, globaldefs.h, checkdefs.F, def_mod.F, extract_obs.F, initial.F, inp_par.F, mod_fourdvar.F, mod_ncparam.F, mod_scalars.F, obs_read.F, obs_write.F, and output.F. Mostly all the changes involve adding VERIFICATION to few CPP directives.
• Corrected a couple of adjoint bugs in ad_set_vbc.F and ini_adjust.F. These bugs only affects 2D applications and 2D optimal perturbations. Many thanks to Andy Moore for finding and correctiong these bugs.
• Renamed routine stats_4dvar.F to a more appropriate stats_modobs.F. This implied changes in is4dvar_ocean.h, s4dvar_ocean.h, w4dpsas_ocean.h, and w4dvar_ocean.h.
• Corrected parallel bugs in grid_coords.F and interpolate.F when computing the initial locations of stations and floats from (lon,lat) data. This bug was reported in the ROMS forum by Wen Long.
• Corrected an array out-of-bounds in nf_fread.F in distributed-memory applications. The scratch variable wrk need to be allocated as Npts+2 since the minimum and maximum values are appended to the end of wrk in routine mp_scatter to reduce additional parallel reduction. Many thanks to John Warner for bringing this to my attention. This one was a hard one.
• Corrected bugs in ecosim.h to fix how pigment biomass is computed. Many thanks to Bronwyn for fixing this one.
• Corrected a typo in wrt_station.F. Many thanks to Jacopo for finding this one.
• Corrected computation of volume-averaged kinetic energy in diag.F, ad_diag.F, rp_diag.F, and tl_diag.F as reported by Bin Zhang and Sasha in the forum.
• Correct several parallel bugs when processing point data in set_data.F, ad_set_data.F, rp_set_data.F, and tl_set_data.F. Many thanks to Rob Hetland for finding this one.

I am considering this blog as a location for my notes on major model advancements and as an additional means to track what has been incorporated into the releases.

Subversion version control.
The idea is that we start with a distributed version of ROMS, roms_v2.2.
This was our starting point and we entered roms_v2.2 into SVN version control. As we make modifications for the sediment advancements, we will make new (internal) releases and use the second dot (which is the third number), such as
roms_sed_v2.2.1, v2.2.2, v2.2.3, etc etc. As more features become available, and as we gain confidence in the third number, we will merge back to the main roms trunk and make roms_v2.3.

*******************************************************************

Starting with roms_v2.2 we have added capabilites for:
– ssw_bbl : bottom boundary layer model to account for increased bottom stress due to wave+currents

– suspended sediment transport : minor additions to account for multiple sediment classes and new input file structures.

– bedload transport : included 2 methods 1) Meyer-Peter Mueller and 2) Soulsby +Damgaard methodologies. Both use bottom stress for waves + currents but the SD method is more applicable to wave environemnts. The MPM method was developed just for currents.

– bed model : 3D bed strucutre that tracks stratigraphy, multiple sed class distributions, resuspension + deposition interaction with suspended sed.

– morphology : as bed evolves, the changing sea floor elevation is used as boundary condition to omega. Entire model knows about changing bed elevation (if activated).

– wetting/drying : this is based on cell face blocking methodology. The user specifies a ‘critical depth.’ During computation, if the total water depth (h+zeta) at a rho point is less than the critical depth, then transport is prevented out of that cell. That is it. Water is always allowed to flow into a cell. Cells with an initial rho_mask = 0 will always be dry.

–surface tke fluxes : there are 2 methods to account for flux of turbulent kinetic energy from breaking surface waves 1) Craig/Banner with Charnok coefficient that is based on a surface wind stress and 2) wave energy dissipation that is based on wave breaking. Both methods are incorporated into the GLS turbulence closure model.

– wave/current interaction : this is currently based on Mellor 2003; 2005 JPO and adds ‘radiation stress terms’ to the momentum equations. These terms add forcing due to momentum flux from waves. The formulation requires information of wave height, length, and direction.

– model coupling : we use the Model Coupling Toolkit to couple ROMS to SWAN (v40.41AB). Users need to install and compile MCT to make the librarires for linking. The model is driven with the ROMS/External/coupling*.in file.

These features now comprise our version roms_sed_v2.2.2.

*******************************************************************
roms_sed_v2.2.3
– minor updates to several test cases

*******************************************************************
roms_sed_v2.2.4
– update to new directory structure of

— src
—-Compilers
—-Lib
—-Master
—-ROMS
——Adjoint
——Bin
——Drivers
——External
——Include
——Modules
——Nonlinear
——Programs
——Representer
——SeaIce
——Tangent
——Utility
—-SWAN

This is code from checkvars.F:

   +84  #ifndef ANA_INITIAL
   +85        get_var(idFsur)=.TRUE.
   +86        get_var(idUbar)=.TRUE.
   +87        get_var(idVbar)=.TRUE.
   +88  # ifdef SOLVE3D
   +89        get_var(idUvel)=.TRUE.
   +90        get_var(idVvel)=.TRUE.
   +91        DO itrc=1,NAT
   +92          get_var(idTvar(itrc))=.TRUE.
   +93        END DO
   +94  # endif
   +95  #endif

I am having trouble in the restart where it reads the time from the restart file, but nothing else. I have a case with #define ANA_INITIAL because I want zero initial flow and zeta (2-D tides), but want to spin up the tides, then restart and save more often. I can’t be the only one who ever wants to restart a test problem. If we ever tackle the problems in the restart, a test problem is an ideal situation to be investigating. I’ll get back to you with a fix…

Problems in restart you ask? Yes, I believe there are some, one being the case of changing history files every say ten saves, then restarting after five or fifteen saves. Try it, you’ll like it.

Edit: I’ve found a fix, hope it works.

A few more mods:

• src/ROMS/Utility/ran_state.F & ran1.F: where necessary, replace integers of default kind with integers of kind i4b (4-byte). This is necessary on Cray T3E, and maybe other platforms, where the default integer is not 4 bytes long.
• src/ROMS/Utility/metrics.F: initialise my_grdmax.

File attached, with any luck