Generalized Stability Theory Algorithms Revisited

[arango]

All the Generalized Stability Theory (GST), adjoint-based algorithms were updated:

    afte_ocean.h      Adjoint Finite Time Eigenmodes (AFT_EIGENMODES)
    fte_ocean.h       Finite Time Eigenmodes (FT_EIGENMODES)
    fsv_ocean.h       Forcing Singular Vectors (FORCING_SV)
    op_ocean.h        Optimal Perturbations (OPT_PERTURBATION)
    so_semi_ocean.h   Stochastic Optimals, Seminorm (SO_SEMI)

• Corrected problem with the convergence of the non-symmetric algorithms that used PARPACK routine pdnaup2.f. It turns out that I introduced a bug when I implemented the check-pointing (restart) procedure. I basically screw-up the ido parameter in this routine. At line 282, we need to have:

            if (ido .eq. 0) then
  
  instead of
  
            if (ido .eq. -2) then
  

  The ido=-2 is used for restart. Without this correction, all the non-symmetric algorithms coverged at the NCV Arnoldi iterations. It took me almost a month to find this bug ...

• The check-pointing CPP option, CHECKPOINTING, is not longer an internal flag. It is removed from globaldefs.h and the User needs to activate it explicitly in order to create/use the GST restart NetCDF file. Recall the the restart switch LrstGST is activated in standard input ocean.in file.

• Corrected a parallel bug when writing the GST restart NetCDF file. The code was hanging-up in a collective MPI communication. The subroutine mp_ncwrite in distributed.F is replaced by mp_ncwrite1d and mp_ncwrite2d to write out 1D/2D state arrays.

• I removed all the calls to the BLAS library routines daxpy and dnrm2, which is used to computed the Eucledean norm in the GST algorithms. The BLAS library only works in parallel with the PGI compiler. We were getting the wrong norms with the ifort and gfortran compilers. This is due to the legacy code style in this library that dimensions arrays as x(1) instead of x(n). The MPI reduction operations were giving incorrect values. In addition, these routines and others have deprecated arithmetic GO TO and were getting a lot of warnings during compilation. Some of these legacy libraries need to be updated for future compilers that do not support old f77 dimension style. This stuff is scary and spend also a lot of time trying to figure out what was going on. To avoid problems in the future, I wrote my own Eucledean norm routines for real vectors (r_norm2) and complex vectors (c_norm2).

• Fixed a bug in ad_ini_fields.F when computing ouput arrays ad_ubar_sol and ad_vbar_sol. We were getting a zero value for ad_ubar and ad_vbar. This was affecting only the GST and adjoint sensitivity algorithms and not the 4D-Var drivers. We were loading data to these arrays after their values have been zero out in the adjoint vertical integration of 3D momentum. We just needed to fill these arrays before of the vertical integration.

• The output Ritz vectors from the GST algorithms can be written in separate NetCDF when the new standard input logical switch LmultiGST is activated:

  ! GST output and check pointing restart parameters.
  
     LmultiGST =  F                               ! one eigenvector per file
       LrstGST =  F                               ! GST restart switch
    MaxIterGST =  500                             ! maximun number of iterations
          NGST =  10                              ! check-pointing interval
  

  The eigenvectors for the AFT_EIGENMODES or FT_EIGENMODES are usually complex. In this case, both the real and imaginary eigenvectors are stored in the same file when LmultiGST is activated. Notice that all the eigenvectors have an output record for the initial or final perturbation, except the eigenvector that have only the adjoint model as a propagator (AFT_EIGENMODES and SO_SEMI). This is a function of the running time window.

• To facilitate the plotting of these eigenvectors from a single or multiple output NetCDF files, several scripts are provided in the ROMS plotting package: plt_afte.sh, plt_fsv.sh, plt_fte.sh, plt_op.sh, and plt_so_semi.sh (see ​src:ticket:493).

• The Forcing Singular Vectors and Stochastic Optimals Vectors are expanded to include the full state (zeta, u, v, temp, salt, ...) and/or surface forcing (sustr, svstr, shflx, ssflux, ...). Several new logical switches were added to standard input ocean.in to set the desired state vector:

  ! Logical switches (TRUE/FALSE) to specify the state variables for
  ! which Forcing Singular Vectors or Stochastic Optimals is required.
  
  Fstate(isFsur) == F                        ! free-surface
  Fstate(isUbar) == F                        ! 2D U-momentum
  Fstate(isVbar) == F                        ! 2D V-momentum
  Fstate(isUvel) == F                        ! 3D U-momentum
  Fstate(isVvel) == F                        ! 3D V-momentum
  Fstate(isTvar) == F F                      ! NT tracers
  
  Fstate(isUstr) == T                        ! surface U-stress
  Fstate(isVstr) == T                        ! surface V-stress
  Fstate(isTsur) == F F                      ! NT surface tracers flux
  

---

WARNINGS:

• Changes were made to the ARPACK/PARPACK libraries!!! The User needs to recompile these libraries to correct the problem with pdnaup2.f described above.

• All the standard input ocean_*.in were changed to add new parameters. This is only relevant to the GST algorithms.

• There is some problems when running these drivers with ifort in the OS Mac environment. This is due to the small stack size in the Mac's. We usually get the following error:

  forrtl: severe (174): SIGSEGV, segmentation fault occurred
  

  This problems will disappear in when I release Phase I of the nesting update in Version 3.5 in few days. Stay tuned, a major ROMS update is coming in your way soon...

---

Personally, this is one of the most difficult updates to ROMS that I have done in years. I have been struggling with this for the last three months. A lot of frustrations with compilers and debuggers. I came across a challenging technical development that I needed to solve in the new nesting infra-structure. I have been thinking about this problem for more than one year.