Custom Query (986 matches)

Corrected the lower time value (Tstr) in get_cycle.F for adjoint applications. Now we have the following statements around line 161:

        tstart=Tmin
        IF (model.eq.iADM) THEN          ! backward time logic
          DO i=2,ntime
            IF ((tstart.lt.mday).and.(mday.le.Tval(i))) THEN
              Tindex=i-1
              Tstr=tstart
              EXIT
            END IF
            tstart=Tval(i)
          END DO
        ELSE                             ! forward  time logic
          DO i=2,ntime
            IF ((tstart.le.mday).and.(mday.lt.Tval(i))) THEN
              Tindex=i-1
              Tstr=tstart
              EXIT
            END IF
            tstart=Tval(i)
          END DO
        END IF

We need to have a different inequality in the IF-statement for backward time logic in adjoint applications. This was affecting the GST drivers when selecting the initial snapshot times to interpolate the forward trajectory. It is weird that we didn't noticed this problem in the 4D-Var algorithms. The reason for that is that the time window between the NLM and ADM is the same. In the GST algorithms, we can have smaller time-window from the forward trajectory and the bug appears...

I also initialized the local LwrtGST switch in the afte_ocean.h driver to allow check-pointing.

There is a parallel issue with the GST drivers. There is an MPI hang-on during the iterations. I am hunting for this elusive parallel (MPI) problem.

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.

A new vertical stretching function, Vstretching = 4, was added and it is denoted as Shchepetkin (2010). This is a major improvement to his 2005 function, Vstretching = 2. The new function, C(s), has a double stretching transformation:

   C(s) = [1 - COSH(theta_s * s)] / [COSH(theta_s) - 1]

with bottom refinement:

   C(s) = [EXP(theta_b * C(s)) - 1] / [1 - EXP(-theta_b)]

Notice that the bottom function is the second stretching of an already stretched transform. The resulting function is continuous with respect the control parameters theta_s and theta_b with a meaningful range of:

   0 <  theta_s <= 10.0
   0 <= theta_b <=  4.0

Due to its functionality and properties Vtransform = 2 and Vstretching = 4 are now the default values for ROMS.

Notice that with Vtransform = 2 the true sigma-coordinate system is recovered as hc goes to infinity. This is useful when configuring applications with flat bathymetry and uniform level thickness. Practically, you can achieve this by setting in ocean.in:

   THETA_S = 0.0d0
   THETA_B = 0.0d0
   TCLINE  = 1.0d+17       (a large number)

Please check ​WikiROMS for updated documentation and plots. For example, in the coarse resolution, North Atlantic application (DAMEE_4) we get:

WARNING:

The new vertical stretching transformation Vtransform = 2 and Vstretching = 4 may have different time-step limitations due to the thickness of the discritized vertical levels. Users need to pay attention to extreme r-factor (rx1) values near the bottom

The following testcases were changed to have the default transformation Vtransform = 2 and Vstretching = 4:

                       OLD (1,1)                      NEW (2,4) 
               theta_s  theta_b  Tcline      theta_s   theta_b  Tcline
              ---------------------------------------------------------

BASIN            3.3      0.0     1400         3.3       0.0     1400
BENCHMARK        5.0      0.4      200         0.0       0.0      400                      
BL_TEST          3.2      0.2        5         5.0       1.5       20
CANYON2D        1d-4      0.0       90         0.0       0.0    1d+16
DOUBLE_GYRE      1.0      1.0       50         0.0       0.0    1d+16
FLT_TEST2D       3.0      0.4      100         0.0       0.0    1d+16
FLT_TEST3D       3.0      0.4      100         0.0       0.0    1d+16
GRAV_ADJ        1d-4      0.0       50         0.0       0.0    1d+16
KELVIN           3.0      0.0       50         0.0       0.0    1d+16
LMD_TEST        1d-4      0.0       50         0.0       0.0    1d+16
RIVERPLUME1      3.0      0.4       50         3.0       0.0       30
RIVERPLUME2      3.0      0.4       50         3.0       0.0       30
SEAMOUNT         5.0      0.4       50         6.5       2.0      100
SOLITON         1d-4      0.0       50         0.0       0.0    1d+16
UPWELLING        3.0      0.0       50         3.0       0.0       25
WINDBASIN       1d-4      0.0       50         1.0       0.0       50

Many thanks to Sasha Shchepetkin for designing this nice and continuous transformation.