;; 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= n_elements(star(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(0)
AMP2= AMP(1)

PHASE1= PHASE(0)
PHASE2= PHASE(1)

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" )

!p.multi= [0,2,3]

;; 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=replicate(1,13)
hplot, PHASE2, fill=1, demarc=[3], levs=indgen(7), c_labels=replicate(1,7)

end