For ocean tidal predictions, once an agreement is reached on the essential lunisolar terms, then the secondorder terms are refined. Early in the last century Doodson catalogued most of these terms:
“Since the midtwentieth century further analysis has generated many more terms than Doodson’s 388. About 62 constituents are of sufficient size to be considered for possible use in marine tide prediction, but sometimes many fewer can predict tides to useful accuracy.”
That’s possibly the stage we have reached in the ENSO model. There are two primary terms for lunar forcing (the Draconic and Anomalistic) cycles, that when mixed with the annual and biannual cycles, will reproduce the essential ENSO behavior. The secondorder effects are the modulation of these two lunar cycles with the Tropical/Synodic cycle. This is most apparent in the modification of the Anomalistic cycle. Although not as important as in the calculation of the Total Solar Eclipse times, the perturbation is critical to validating the ENSO model and to eventually using it to make predictions.
The variation in the Anomalistic period is described at the NASA Goddard eclipse page. They provide two views of the variation, a timedomain view and a histogram view.


Time domain view 
Histogram view 
Since NASA Goddard doesn’t provide an analytical form for this variation, we can see if the ENSO Model solver can effectively match it via a bestfit search to the ENSO data. This is truly an indirect method.
First we start with a parametric approximation to the variation, described by a pair of successive frequency modulated (and fullwave rectified) terms that incorporate the Tropicalmodified term, wm. The Anomalistic term is wa.
COS(wa*t+pa+c_1*ABS(SIN(wm*t+k_1*ABS(SIN(wm*t+k_2))+c_2)))
This can generate the cusped behavior observed, but the terms pa, c_1, c_2, k_1, and k_2 need to be adjusted to align to the NASA model. The solver will try to do this indirectly by fitting to the 18801950 ENSO interval.
Plotting in RED the Anomalistic time series and the histogram of frequencies embedded in the ENSO waveform, we get:


Time domain view of model 
Histogram view of model 
This captures the histogram view quite well, and the timedomain view roughly (in other cases it gives a better cusped fit). The histogram view is arguably more important as it describes the frequency variation over a much wider interval than the 3year interval shown.
What would be even more effective is to find the correct analytical representation of the Anomalistic frequency variation and then plug that directly into the ENSO model. That would provide another constraint to the solver, as it wouldn’t need to spend time optimizing for a known effect.
Yet as a validation step, the fact that the solver detects the shape required to match the variation is remarkable in itself. The solver is obviously searching for the forcing needed to produce the ENSO waveform observed, and happens to use the precise parameters that also describe the secondorder Anomalistic behavior. That could happen by accident but in that case there have been too many happy accidents already, i.e. period match, LOD match, Eclipse match, QBO match, etc.