// Example HOPS script #7
// harmonic: Harmonic analysis (amplitude/phase) of the annual cycle 

// Open 1 year of prognostic ocean data (24 samples at half-month intervals)
// for a 3-degree NCAR Ocean Model integration
// (converted to CSM netCDF format using Frank Bryan"s filters)

hopen, "/data/large1/CSM/3x3/G03-4.24/Y0989_0990_p.nc";

// Get surface (z=0.) temperature (T) 
SST= hget("T",z=0.);

// Compute standard deviation
SDEV= hsub( hop(SST, "-", hsub(SST,t="avg")), t="rms" );
SDEV.name= "Tsdev";

// Number of time samples
nt= numberof(*(SST.time));

// No. of days per year
ndpyr= hattr(SST, "time:days_per_year" );

// Compute FFT
FFT_SST= hfft(SST,"t",nocheck=1);

// Extract mean from FFT output
MEAN= hop( "real", hop( hsub( FFT_SST, t=0.), "/", sqrt(float(nt)) ) );
MEAN.name= "Tmean";

// Extract annual, semiannual harmonics from FFT output
HARMONIC= hop( hsub(FFT_SST, t=[1./ndpyr,2./ndpyr]), "/", 0.5*sqrt(float(nt)) );

// Compute amplitude
AMP= hop( "abs", HARMONIC );

// Compute phase of annual cycle (in radians, -pi to pi)
PHASE= hop("-", hop( hop("imaginary",HARMONIC), "atan", 
                     hop("real",HARMONIC) ));

// Ensure that phase is positive
PHASE= hop( PHASE, "+", hop( hop( PHASE,"<",0.), "*", 2*pi()) );

// Split annual and semi-annual harmonics into separate slabs
AMP= hsplit(AMP,"t");
PHASE= hsplit(PHASE,"t");

AMP1= AMP(1);
AMP2= AMP(2);

PHASE1= PHASE(1);
PHASE2= PHASE(2);

PHASE1.name= "Tphase1";
PHASE2.name= "Tphase2";

// Convert phase to month when maximum SST occurs
PHASE1= hop( PHASE1, "*", 6./pi(), units="months" );
PHASE2= hop( PHASE2, "*", 3./pi(), units="months" );

// Plot mean
hplot, MEAN, fill=1, levs="";

// Plot standard deviation
hplot, SDEV, fill=1, levs="";

// Plot amplitudes
hplot, AMP1, fill=1, levs="";
hplot, AMP2, fill=1, levs="";

// Plot phases
hplot, PHASE1, fill=1, demarc=[6], levs=indgen(13), c_labels=array(1,13);
hplot, PHASE2, fill=1, demarc=[3], levs=indgen(7), c_labels=array(1,7);