Generalized Stability Theory Algorithms Revisited

[arango]

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

    afte_ocean.h      Adjoint Finite Time Eigenmodes (AFTE_EIGENMODES)
    fte_ocean.h       Finite Time Eigenmodes (FTE_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 are flipped so the the most significant is written first instead of last (unlimited dimension) in the NetCDF files.

• 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.