ROMS/TOMS Developers

I introduced the following updates to the code:

• Add new time-averaged quantities which are needed to compute eddy fluxes. For example, the tracer eddy flux is:

  <H u T / n> = <H u / n> <T> + <H’ u’ t’ / n>
  <H v T / m> = <H v / m> <T> + <H’ v’ t’ / m>

  we are interested in the last term which can be computed as:

  <H’ u’ t’ / n> = <H u / n> – <H u T / n> <T>
  <H’ v’ t’ / m> = <H v / m> – <H v T / m> <T>

  where the quantities inside <…> indicate time-averaged. In model variables, the tracer eddy flux is:

  <H’ u’ t’ / n> = <Huon*t> – <Huon> <t>
  <H’ v’ t’ / m> = <Hvom*t> – <Hvom> <t>

  all the three time-averaged quantities in the right-hand-size are now saved in the averages NetCDF file and can be activated with cpp option AVERAGES_QUADRATIC. Many thanks to John Wilkin for helping us to implement this option.

• Modified cppdefs.h in such a way that each application has it own cpp definitions include file:

  /*
  **-----------------------------------------------------------------------------
  **  Include the appropriate header file for IDEALIZED APPLICATION.
  **-----------------------------------------------------------------------------
  */
  
  #if defined A4DVAR_TOY
  # include "a4dvar_toy.h"
  #elif defined BASIN
  # include "basin.h"
  #elif defined BENCHMARK1 || defined BENCHMARK2 || defined BENCHMARK3
  # include "benchmark.h"
  #elif defined  BIO_TOY
  # include "bio_toy.h"
  #elif defined BL_TEST
  # include "bl_test.h"
  #elif defined CANYON_A || defined CANYON_B
  # include "canyon.h"
  #elif defined CHANNEL_NECK
  # include "channel_neck.h"
  #elif defined COUPLING_TEST
  # include "coupling_test.h"
  #elif defined DOUBLE_GYRE
  # include "double_gyre.h"
  #elif defined ESTUARY_TEST
  # include "estuary_test.h"
  #elif defined FLT_TEST
  # include "flt_test.h"
  #elif defined GRAV_ADJ
  # include "grav_adj.h"
  #elif defined INNER_PRODUCT
  # include "inner_product.h"
  #elif defined KELVIN
  # include "kelvin.h"
  #elif defined LAB_CANYON
  # include "lab_canyon.h"
  #elif defined LAKE_SIGNELL
  # include "lake_signell.h"
  #elif defined LMD_TEST
  # include "lmd_test.h"
  #elif defined OVERFLOW
  # include "overflow.h"
  #elif defined RIVERPLUME1 || defined RIVERPLUME2
  # include "riverplume.h"
  #elif defined SEAMOUNT
  # include "seamount.h"
  #elif defined SED_TEST1
  # include "sed_test.h"
  #elif defined SED_TOY
  # include "sed_toy.h"
  #elif defined SOLITON
  # include "soliton.h"
  #elif defined UPWELLING
  # include "upwelling.h"
  #elif defined WEDDELL
  # include "weddell.h"
  #elif defined WINDBASIN
  # include "windbasin.h"
  
  /*
  **-----------------------------------------------------------------------------
  **  Include the appropriate header file for REALISTIC APPLICATION.
  **-----------------------------------------------------------------------------
  */
  
  #elif defined ADRIA02
  # include "adriatic.h"
  #elif defined CBLAST
  # include "cblast.h"
  #elif defined DAMEE_4
  # include "damee.h"
  #elif defined EAC_4 || defined EAC_8
  # include "eac.h"
  #elif defined IAS
  # include "ias.h"
  #elif defined NATL
  # include "natl.h"
  #elif defined NENA
  # include "nena.h"
  #elif defined NJ_BIGHT
  # include "nj_bight.h"
  #elif defined NPACIFIC
  # include "npacific.h"
  #elif defined SCB
  # include "scb.h"
  #elif defined SW06_COARSE || defined SW06_FINE
  # include "sw06.h"
  #elif defined USWEST
  # include "uswest.h"
  #else
    CPPDEFS - Choose an appropriate ROMS application.
  #endif
  
    

  This new design is cleaner since cppdefs.h was getting very big. I have noticed that users like to put a lot of comments in their file to document their application. The user needs now to activate the appropriate application option in cppdefs.h (for example, #define UPWELLING) and make the appropriate changes to its respective definitions file (upwelling.h). All the include files are located in the Include sub-directory. The user needs to create an include file for each different application. I wanted to do this for very long time. The feedback that I got from several of you about this was very positive.

• Several of you have been requesting for a perfect restart option in ROMS. I started coding this option. A perfect restart option is provided for the basic kernel. It can be activated with cpp option PERFECT_RESTART. This option is still under testing. I have tested this option with the upwelling example. It is run for 20 time-steps. The perfect restart fields are written at time-step 10. The RMS errors at time-step 20 between restarted and control run are:

  
                            RMS               RMS               RMS
     zeta              2.73101330e-06    2.73138360e-06    2.73142573e-06
     ubar              1.50954195e-06    1.50982080e-06    1.50982771e-06
     vbar              1.47805759e-07    1.47870988e-07    1.47942021e-07
     temp              1.12177271e-05    1.30822796e-05    1.40771552e-05
     salt              2.70888331e-05    2.90445209e-05    3.22833452e-05
     u                 1.50989621e-06    1.54476679e-06    1.54760557e-06
     v                 5.84640682e-07    1.13022193e-06    1.17730862e-06
                         ANA_VMIX          GLS_MIXING        MY25_MIXING

  However, if PERFECT_RESTART is not activated, we get:

  
                            RMS               RMS               RMS
     zeta              1.71578595e-09    5.68765614e-09    2.07986550e-09
     ubar              3.91111746e-09    4.38520664e-09    4.21003794e-09
     vbar              9.83973064e-11    2.34650383e-09    1.13328922e-10
     temp              4.01893636e-07    7.00800631e-05    5.55409335e-07
     salt              0.00000000E+00    0.00000000E+00    0.00000000E+00
     u                 9.59030197e-08    1.05812267e-06    2.71992162e-07
     v                 2.89748763e-07    7.39400620e-07    8.01745441e-07
                         ANA_VMIX          GLS_MIXING        MY25_MIXING

  Obviously, this option is not working yet. Two new routines were introduced to achieve the perfect restart in a generic way: nf_read4d.F and nf_write4d.F.

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.