HRET

# High Resolution Empirical Tides

# Computing tidal predictions with the HRET models

The HRET models are described [here][http://web.cecs.pdx.edu/~zaron/pub/HRET.html].
The NetCDF files contain real and imaginary parts of the harmonic constants. Thus,
the predicted tide at a particular time, `t`, may be computed from `h(t) = re*cos(omega*t) + im*sin(omega*t)`,
where time is provided as a real-valued Julian Date, and the frequency of the tide is, `omega`.
For more precision, the nodal modulations of amplitude and phase may be included, and `omega*t` may
be computed directly from the mean longitudes of the relevant astronomical quantities.

The software in this repository contains different versions of this computation, which
should allow you to compute tidal predictions of individual tidal components, e.g., M2, or from
the sum of all the components in the HRET model.

# Sub-directories:

 - ./hret_tide/ -- Remko Scharoo's fortran software for the Jason-CS Preprocessor. Received 2019-03-08.
 This sofware contains everything you need to compute the arguments of the cosine and sine functions
 for the given times.

 - ./hret_vel/ -- My julia-language software for computing tidal currents at the ocean surface. This software
 calls my harmonic analysis library `jd_funcs.jl`, in the [HA repository][https://hachi.cee.pdx.edu/fossil/HA].

# Tidal periods

| Darwin Symbol | Period (days) |
---------------------------------
| M2            | 0.5175250544  |
| S2            | 0.5000000000  |
| K1            | 0.9972695425  |
| O1            | 1.0758058718  |
| MA2           | 0.5182593837  |
| MB2           | 0.5167927854  |
---------------------------------

 - While I have made every effort to provide accurate information, I urge you to double check the numerical
results you compute with this software. Please contact me if you find mistakes in the information presented here.

 - Be aware that the HRET models are models of the *internal tide*. The predicted elevations are typically only
 a few centimeters. For ocean tide predictions at ports and harbors, please
 see [the NOAA website][http://tidesandcurrents.noaa.gov/].

########################################################
# Makefile for the HRET internal tide model source
########################################################
# (c) Remko Scharroo - EUMETSAT

#-----------------------------------------------------------------------
# Additional tools
#-----------------------------------------------------------------------
LN_S = ln -s -f
RANLIB = ranlib
FC = $(shell nc-config --fc)
FFLAGS := $(shell nc-config --fflags)
LDLIBS := $(shell nc-config --flibs)

#-----------------------------------------------------------------------
# Special rules
#-----------------------------------------------------------------------
# .f90 rules
%.o:	%.f90
	$(COMPILE.f) $(OUTPUT_OPTION) $<
%:	%.f90
	$(LINK.f) $^ $(LOADLIBES) $(LDLIBS) -o $@

# Cancel default .mod rules
%:	%.mod
%.o:	%.mod

# Create new .mod rule
# Do not add %.f90 dependency because %.mod retain timestamp when not altered
%.mod:
	$(COMPILE.f) -o $*.o $*.f90

do_hrettidetest:	hrettidetest
	hrettidetest HRET_v8.1_compressed.nc

hrettidetest:	tides_hret.o tides_aux.o

tides_hret.o:	tides_aux.mod

clean:
	rm -f *.o *.mod hrettidetest

tar:	clean
	(cd ..; tar -cvzhf hret_tide.tar.gz hret_tide)

program rewrite_hret

! Rewrite the HRET internal tide grid file such that:
! - It uses netCDF-4 internal compression
! - The maps are expanded to 90 latitude (from 66)
! - Values are stored as short instead of double
! - Some additional attributes are added conforming to CF-1.7
!-
! 19-Oct-2018 - Created by Remko Scharroo (EUMETSAT)
!-----------------------------------------------------------------------
use typesizes
use netcdf
use rads_netcdf
use rads_misc

! Initialisation

integer(fourbyteint), parameter :: nx=7200, ny=3601, deflate=4
real(eightbytereal), parameter :: x0=0d0, x1=359.95d0, y0=-90d0, y1=90d0, dx=0.05d0, dy=0.05d0, scale_factor=1d-4
real(eightbytereal) :: x(nx), y(ny)
character(80) :: filenm1, filenm2, units, long_name, varnm
integer(fourbyteint) :: ncid1, ncid2, dimid(2), id, i, ny0, ny1, kx, ky, nvar, natt
integer(twobyteint) :: fill_value=-32767-1
real(eightbytereal), allocatable :: orig(:,:)
integer(twobyteint), allocatable :: new(:,:)
integer(fourbyteint), allocatable :: varid(:)

! Read input and output filename from command line

call getarg(1,filenm1)
call getarg(2,filenm2)

! Load the original grid

call nfs(nf90_open(filenm1,nf90_nowrite,ncid1))
call nfs(nf90_inquire(ncid1,nvariables=nvar,nattributes=natt))
allocate (varid(nvar))
call nfs(nf90_inq_dimid(ncid1,'latitude',dimid(1)))
call nfs(nf90_inq_dimid(ncid1,'longitude',dimid(2)))
call nfs(nf90_inquire_dimension(ncid1,dimid(1),len=ky))
call nfs(nf90_inquire_dimension(ncid1,dimid(2),len=kx))
call nfs(nf90_inq_varid(ncid1,'latitude',varid(1)))
call nfs(nf90_inq_varid(ncid1,'longitude',varid(2)))
call nfs(nf90_get_var(ncid1,varid(1),y(:ky)))
call nfs(nf90_get_var(ncid1,varid(2),x(:kx)))

! Allocate memory

allocate (orig(kx,ky),new(nx,ny))
ny0 = nint((y(1)-y0)/dy+1)
ny1 = nint((y(ky)-y0)/dy+1)

! Open the new grid

call nfs(nf90_create(filenm2,nf90_write+nf90_nofill+nf90_netcdf4,ncid2))

! Copy variables and attributes

call nfs(nf90_put_att(ncid2,nf90_global,'Conventions','CF-1.7'))
call nf90_def_axis(ncid2,'longitude','longitude','degrees_east',nx,x0,x1,dimid(1),varid(1), &
	'X','longitude',deflate)
call nf90_def_axis(ncid2,'latitude','latitude','degrees_north',ny,y0,y1,dimid(2),varid(2), &
	'Y','latitude',deflate)
i = 2
do id = 1,nvar
	call nfs(nf90_inquire_variable(ncid1,id,name=varnm))
	if (varnm == 'latitude' .or. varnm == 'longitude') cycle
	i = i + 1
	call nfs(nf90_def_var(ncid2,varnm,nf90_int2,dimid,varid(i)))
	call nfs(nf90_get_att(ncid1,id      ,'longname',long_name))
	call nfs(nf90_put_att(ncid2,varid(i),'long_name',long_name))
	call nfs(nf90_get_att(ncid1,id      ,'units',units))
	call nfs(nf90_put_att(ncid2,varid(i),'units',units))
	call nfs(nf90_put_att(ncid2,varid(i),'scale_factor',scale_factor))
	call nfs(nf90_put_att(ncid2,varid(i),'_FillValue',fill_value))
	call nfs(nf90_def_var_deflate(ncid2,varid(i),1,deflate,deflate))
	write (*,*) id,i,varid(i),varnm
enddo
do id = 1,natt
	call nfs(nf90_inq_attname(ncid1,nf90_global,id,varnm))
	call nfs(nf90_copy_att(ncid1,nf90_global,varnm,ncid2,nf90_global))
enddo
call nfs(nf90_put_att(ncid2,nf90_global,'compression','compress_hret '//trim(filenm1)//' '//trim(filenm2)))
call nfs(nf90_enddef(ncid2))

! Write coordinate axes

call nf90_put_axis(ncid2,varid(1))
call nf90_put_axis(ncid2,varid(2))

! Write grids

new = 0
i = 2
do id = 1,nvar
	call nfs(nf90_inquire_variable(ncid1,id,name=varnm))
	if (varnm == 'latitude' .or. varnm == 'longitude') cycle
	i = i + 1
	call nfs(nf90_get_var(ncid1,id,orig))
	new(:,ny0:ny1) = nint2(orig/scale_factor)
	write (*,*) id,i,varid(i),varnm
	call nfs(nf90_put_var(ncid2,varid(i),new))
enddo
call nfs(nf90_close(ncid1))
call nfs(nf90_close(ncid2))

deallocate (orig,new,varid)

end program rewrite_hret

program hrettidetest
!
! This program tests the hrettide routine based on
! information provided by Ed Zaron
!
! Syntax: hrettidetest <pathname>
!
! where <pathname> is the pathname of the HRET NetCDF file
!-
! $Log: hrettidetest.f90,v $
! Revision 1.2  2019/03/08 09:57:29  rads
! - Introduced internal tide model HRET
!
! (c) Remko Scharroo - EUMETSAT
!-----------------------------------------------------------------------
use tides_hret
use typesizes

! Initialise

integer, parameter :: nt = 3
integer :: i
logical :: passed
character(len=512) :: pathname
type(hrettideinfo) :: info
type :: ref
	real(eightbytereal) :: time, lat, lon, tide_comp(6)
endtype
real(eightbytereal) :: time, tide, tide_comp(6)
type(ref) :: t(nt)

! Test cases

t(1) = ref(1046250908.335094d0, 18.948867d0, 149.491272d0, &
	(/-0.004648d0, -0.000328d0, 0.004537d0, 0.001753d0, 0.001196d0, 0.000073d0/))
t(2) = ref(1048578767.430392d0, 17.925634d0, 149.995549d0, &
	(/-0.001640d0, 0.001641d0, 0.000251d0, -0.000959d0, 0.000565d0, -0.001284d0/))
t(3) = ref(1051079248.994491d0, 18.960929d0, 149.837059d0, &
	(/-0.015131d0, -0.000893d0, -0.001623d0, 0.001529d0, 0.001067d0, -0.003253d0/))

! Load the model

call getarg(1, pathname)
call hrettideinit(pathname, info)

! Run the tests

do i = 1, nt
	time = (t(i)%time - 2446066.5d0) * 86400d0
	call hrettide(info, t(i)%time, t(i)%lat, t(i)%lon, tide, tide_comp)
	passed = all(abs(tide_comp - t(i)%tide_comp) < 1d-4)
	write (*,600) i, t(i)
	write (*,601) i, passed, tide_comp
enddo
600 format (i2,f20.6,2f12.6,6f10.6)
601 format (i2,l44,6f10.6)

! Deallocate memory

call hrettidefree(info)

end program hrettidetest

module tides_aux
contains

!*mean_longitudes -- Compute 5 principle astronomical mean longitudes
!+
subroutine mean_longitudes (utc, s, h, p, np, pp)
use typesizes
real(eightbytereal), intent(in) :: utc
real(eightbytereal), intent(out) :: s, h, p, np, pp
!
! This routine computes the five principle mean longtitudes used to compute
! astronomical tides.
!
! Input is the UTC time in days since 1985. This is then converted to
! Julian centuries since 1.5 Jan 2000. Based on the mean longitudes and rates
! at that time, the mean longitudes at the current time are computed.
!
! Output are the mean longitudes of the moon, sun, lunar perigee, ascending
! lunar node, and solar perigee.
!
! Note: This routine uses time in UT and does not distinguish
! between the subtle differences of UTC, UT1, etc. However, this is
! more than accurate enough for our purposes.
! The formalae for the mean longitudes depend on dynamic time (DT).
! This routine assumes DT - UT = 63.48 seconds (2000).
!
! Reference:
! Jean Meeus, Astronomical Algorithms, 2nd ed., 1999.
!
! Arguments:
! utc     : UTC time in days since 1.0 January 1985
! s       : mean longitude of the moon (degrees)
! h       : mean longitude of the sun (degrees)
! p       : mean longitude of the lunar perigee (degrees)
! np      : negative of mean longitude of the ascending lunar node (degrees)
! pp      : mean longitude of the solar perigee (degrees)
!-
! $Log: tides_aux.f90,v $
! Revision 1.3  2019/03/08 09:57:29  rads
! - Introduced internal tide model HRET
!
! Revision 1.2  2017/12/04 15:03:32  rads
! - Update comments only
!
! Revision 1.1  2017/10/12 09:17:24  rads
! - Introduces tides_aux.f90 with generic tide routines, like mean_longitudes
!
! (c) Remko Scharroo - EUMETSAT
!-----------------------------------------------------------------------
real(eightbytereal) :: t

t  = (utc - 5478.4993d0) / 36525d0					! Julian centuries since 1.5 Jan 2000
s  = modulo(218.3164477d0 + 481267.88123421d0 * t, 360d0)	! Mean longitude of moon
h  = modulo(280.4662556d0 +  36000.76983081d0 * t, 360d0)	! Mean longitude of sun
p  = modulo( 83.3532465d0 +   4069.0137287d0  * t, 360d0)	! Mean longitude of lunar perigee
np = modulo(234.95548d0   +   1934.136261d0   * t, 360d0)	! Mean longitude of ascending lunar node
pp =        282.94d0      +      1.7192d0     * t		! Mean longitude of solar perigee

end subroutine mean_longitudes

end module tides_aux

module tides_hret
use typesizes
integer(fourbyteint), parameter, private :: nw = 6
real(eightbytereal), parameter, private :: pi = 4d0*atan(1d0), rad = pi/180d0
character(3), parameter, private :: wave(nw) = (/ 'M2 ', 'K1 ', 'S2 ', 'O1 ', 'MA2', 'MB2' /)
type :: hrettideinfo
	integer(fourbyteint) :: nx, ny, undef(nw)
	integer(twobyteint), allocatable :: p_re(:,:,:), p_im(:,:,:)
	real(eightbytereal) :: xmin, xmax, dx, ymin, ymax, dy, fact(nw), arg(nw), f(nw), u(nw), t_nodal, nan, wmin
endtype

contains

!*hrettide -- Compute internal tides according to HRET model
!+
subroutine hrettide (info, utc, lat, lon, tide, tide_comp)
use tides_aux
type(hrettideinfo), intent(inout) :: info
real(eightbytereal), intent(in) :: utc, lat, lon
real(eightbytereal), intent(out) :: tide
real(eightbytereal), intent(out), optional :: tide_comp(nw)
!
! Compute internal tidal height for given time and location from grids
! of harmonic coefficients of the HRET model. The present version
! uses 6 constituents (M2, K1, S2, O1, MA2, MB2).
! This routine is heavily based on the routine PERTH3 by Richard Ray.
!
! The input files can be found in $ALTIM/data/<name>, where $ALTIM
! is an environment variable and <name> is the argument in the
! call to hrettideinit. The grids are in NetCDF format.
!
! To initialize the computation, the function hrettideinit should be
! called first. It allocates the appropriate amount of memory and
! loads the grids into memory. To release the memory for further
! use, call hrettidefree.
!
! Longitude and latitude are to be specified in degrees; time in UTC
! seconds since 1 Jan 1985. All predicted tides are output in meters.
!
! The computation of nodal periods and mean longitudes is optimized
! for the period 1980-2020.
!
! Input arguments:
!  info     : Structure initialised by hrettideinit
!  utc      : UTC time in seconds since 1 Jan 1985
!  lat      : Latitude (degrees)
!  lon      : Longitude (degrees)
!
! Output arguments:
!  tide     : Predicted internal tide (m)
!  tide_comp: (Optional) Predicted internal tide, separate components (m)
!-
! $Log: hrettide.f90,v $
! Revision 1.2  2019/03/08 09:57:29  rads
! - Introduced internal tide model HRET
!
! (c) Remko Scharroo - EUMETSAT
!-----------------------------------------------------------------------
real(eightbytereal) :: slon
real(eightbytereal), parameter :: mjd85 = 46066d0
integer(fourbyteint) :: istat
complex(eightbytereal) :: val(nw), w(nw)

! When called with invalid time, reset time reference and return

if (.not.(utc < 1d20)) then
	info%t_nodal = 1d30
	tide_comp = info%nan
	tide = info%nan
	return
endif

! If latitude out of range, bail out

if (isnan(lat) .or. isnan(lon) .or. lat < info%ymin .or. lat > info%ymax) then
	tide_comp = info%nan
	tide = info%nan
	return
endif

! Limit longitude to interval of grids

slon = lon
if (slon > info%xmax) then
	slon = slon - 360d0
else if (slon < info%xmin) then
	slon = slon + 360d0
endif

! Compute basic astronomical mean longitudes

call hret_astron (utc)

! Compute all internal tide components

call hret_interp (slon, lat, val, istat)
if (istat == 0) then
	tide_comp = info%nan
	tide = info%nan
else
	w = expj(-(info%arg+info%u)*rad)
	tide_comp = info%f * dble(w*val)
	tide = sum(tide_comp)
endif

contains

!-----------------------------------------------------------------------
! Compute the basic astronomical mean longitudes  s, h, p, N.

subroutine hret_astron (time)
real(eightbytereal), intent(in) :: time
!
! These formulae are for the period 1980 - 2020, and were derived
! from the ASTRO5 code in perth3.f by Richard Ray and are based on:
! Jean Meeus, Astronomical Algorithms, 2nd ed., 1998.
!
! Note 1: This routine uses time in UT and does not distinguish
! between the subtle differences of UTC, UT1, etc. However, this is
! more than accurate enough for our purposes.
! The formalae for the mean longitudes depend on dynamic time (DT).
! This routine assumes DT - UT = 63.48 seconds (2000).
!
! Note 2: The argument w is -n', i.e. w is decreasing with time.
!
! TIME is UTC in seconds since 1985.
!
real(eightbytereal) :: t1,t2,s,h,p,w,pp,sinn,sin2n,cosn,cos2n,t

t = time/86400d0
t1 = modulo (t, 1d0) * 360d0
t2 = 2d0*t1

call mean_longitudes (t, s, h, p, w, pp)
w = -w ! Change n' into w = -n'

info%arg(1) = t2 + 2d0*h - 2d0*s    ! M2
info%arg(2) = t1 + h + 90d0         ! K1
info%arg(3) = t2                    ! S2
info%arg(4) = t1 + h - 2d0*s - 90d0 ! O1
info%arg(5) = info%arg(1) - h       ! MA2
info%arg(6) = info%arg(1) + h       ! MB2

! Determine nodal corrections f and u (only when more than 1 day is passed since last time)

if (abs(time-info%t_nodal) > 86400d0) then
	info%t_nodal = time
	sinn = sin(w*rad)
	cosn = cos(w*rad)
	sin2n = sin(2*w*rad)
	cos2n = cos(2*w*rad)

	info%f( 1) = 1.000d0 - 0.037d0*cosn                 ! M2
	info%f( 2) = 1.006d0 + 0.115d0*cosn - 0.009d0*cos2n ! K1
	info%f( 4) = 1.009d0 + 0.187d0*cosn - 0.015d0*cos2n ! O1

	info%u( 1) = -2.1d0*sinn               ! M2
	info%u( 2) = -8.9d0*sinn + 0.7d0*sin2n ! K1
	info%u( 4) = 10.8d0*sinn - 1.3d0*sin2n ! O1
endif
end subroutine hret_astron

!-----------------------------------------------------------------------
! Interpolates a value from a grid of data at the desired location.
! Interpolation is bilinear.

subroutine hret_interp (x, y, val, istat)
real(eightbytereal), intent(in) :: x, y
complex(eightbytereal), intent(out) :: val(:)
integer(fourbyteint), intent(out) :: istat
integer(fourbyteint) :: kw, i, i0, ii, j, j0, jj
real(eightbytereal) :: wtot, xij, yij, weight(2,2)

istat = 0

! Compute indices for desired position

xij = (x - info%xmin) / info%dx
i0 = floor(xij)
xij = xij - i0

yij = (y - info%ymin) / info%dy
j0 = min(int(yij),info%ny-2) ! Use int() so we do not drop below zero
yij = yij - j0

! Set corner weights

weight(1,1) = (1 - xij) * (1 - yij)
weight(1,2) = (1 - xij) *      yij
weight(2,1) =      xij  * (1-  yij)
weight(2,2) =      xij  *      yij

! Sum up the weighted values for the corners

wtot = 0d0
val = (0d0,0d0)

do i = 1,2
	ii = modulo(i0 + i - 1, info%nx) + 1
	jloop: do j = 1,2
		jj = j0 + j
! Check if the constituents are available:
		if (any(info%p_re(ii,jj,:) == info%undef)) cycle jloop

		istat = istat+1
		wtot = wtot + weight(i,j)

		forall (kw = 1:nw) val(kw) = val(kw) + complex(info%p_re(ii,jj,kw),info%p_im(ii,jj,kw)) * weight(i,j)
	enddo jloop
enddo

if (istat == 0 .or. wtot <= info%wmin) then
	istat = 0
	return
endif
val = val * info%fact / wtot

end subroutine hret_interp

elemental function expj (x)
real(eightbytereal), intent(in) :: x
complex(eightbytereal) :: expj
expj = dcmplx(cos(x),sin(x))
end function expj

end subroutine hrettide

!&hrettideinit -- Initialize GOT tide model
!+
subroutine hrettideinit (pathname, info)
use netcdf
character(*), intent(in) :: pathname
type(hrettideinfo), intent(inout) :: info
!
! Allocate memory for HRET tide modelling and read grids into memory.
!
! Input argument:
!  name     : Pathname of the NetCDF file containing the HRET tide model
!
! Output argument:
!  info     : Structure initialised by hrettideinit
!-
integer(fourbyteint) :: istat, ncid, varid, kw

! Produce log info

write (0,600) trim(pathname)
600 format ('(Loading HRET tide: ',a,')')

! Open and verify netCDF file

if (nf90_open(pathname,nf90_nowrite,ncid) /= nf90_noerr) stop 'hrettideinit: Unable to open file'
if (nf90_inquire_dimension(ncid,1,len=info%nx) + &
    nf90_get_att(ncid,1,'valid_min',info%xmin) + &
    nf90_get_att(ncid,1,'valid_max',info%xmax) /= nf90_noerr) stop 'hrettideinit: Error in longitude dimensions'
if (nf90_inquire_dimension(ncid,2,len=info%ny) + &
    nf90_get_att(ncid,2,'valid_min',info%ymin) + &
    nf90_get_att(ncid,2,'valid_max',info%ymax) /= nf90_noerr) stop 'hrettideinit: Error in latitude dimensions'

! Allocate memory

if (allocated(info%p_re) .or. allocated(info%p_im)) stop 'hrettideinit: Use hrettidefree first'
allocate (info%p_re(info%nx,info%ny,nw), info%p_im(info%nx,info%ny,nw), stat=istat)
if (istat /= 0) stop 'hrettideinit: Not able to allocate memory'

! Initialisations

info%dx = (info%xmax-info%xmin) / (info%nx-1)
info%dy = (info%ymax-info%ymin) / (info%ny-1)
info%fact = 1d0
info%wmin = 0.5d0
info%f = 1d0
info%u = 0d0
info%t_nodal = 1d30
info%nan = 0d0
info%nan = info%nan / info%nan

! Read individual grids

do kw = 1,nw
	if (nf90_inq_varid(ncid,trim(wave(kw))//'re',varid) /= nf90_noerr) stop 'hrettideinit: Missing wave: '//trim(wave(kw))
	if (nf90_get_var(ncid,varid,info%p_re(:,:,kw)) /= nf90_noerr) stop 'hrettideinit: Missing wave: '//trim(wave(kw))
	if (nf90_inq_varid(ncid,trim(wave(kw))//'im',varid) /= nf90_noerr) stop 'hrettideinit: Missing wave: '//trim(wave(kw))
	if (nf90_get_var(ncid,varid,info%p_im(:,:,kw)) /= nf90_noerr) stop 'hrettideinit: Missing wave: '//trim(wave(kw))
	if (nf90_get_att(ncid,varid,'scale_factor',info%fact(kw)) /= nf90_noerr) info%fact(kw) = 1d0
enddo
if (nf90_close(ncid) /= nf90_noerr) return

end subroutine hrettideinit

!&hrettidefree -- Free up space allocated by hrettideinit
!+
subroutine hrettidefree (info)
type(hrettideinfo), intent(inout) :: info
!
! This routine frees up memory allocated by a previous call to hrettideinit
!
! Argument:
!  info : Struct initialised by hrettideinit
!-
if (allocated(info%p_re)) deallocate (info%p_re)
if (allocated(info%p_im)) deallocate (info%p_im)
end subroutine hrettidefree

end module tides_hret

#!/home/ezaron/opt/bin/julia

# 2019-05-20:
#
#  This code is based on namako:NASA-Tides/Julia/driftersx_forProposal.jl.
#
#  Goals:
#    - upgrade to julia-1.*
#    - generate tide current predictions to accompany drifter files.
#

if (length(ARGS) == 0)
    FIG = true
    fdir = "../Drifters/Ver1.02/"
    fu = "driftertrajGPS_1.02.nc"
else
    FIG = false
    fx = ARGS[1]
    m = match(r"(.+)(/.+\.nc)",fx)
    fdir = m[1] * "/"
    fu   = m[2][2:end]
end

println("drifters_predict.jl shall be working on : ",fdir * fu)

include("edzlib.jl")

if (FIG)
    using Winston
    using Colors
end
using NetCDF
using Interpolations
using MAT
using Printf

#include("/home/ezaron/WeddellSea/Julia/edznc.jl")
include("/home/ezaron/HA/iho_func.jl")
using .HAmod

#gamma = parse(Float64,ARGS[1])
gamma = 0.

# To start, just load the hc and see if the currents look okay:
cid = "M2"
#froot = "/home/ezaron/NASA-Tides/HRET_coeus_results/"
#fn    =  froot * cid * "_IZ_global_v7.nc"
cidvec = ["M2" "K1" "S2" "O1" "MA2" "MB2"]
#cidvec = ["M2" "K1" "S2" "O1"]
#cidvec = ["M2" "K1"]
nc = length(cidvec)
froot  = "/home/ezaron/NASA-Tides/Julia/IZ/"
f1     = "HRET_v8.1"
fhret     = froot * "HRET_v8.1" * ".nc"

M      = Array{Dict{String,Any},1}(undef,length(cidvec))

cid = cidvec[1]

# So that the interpolators wrap properly with the periodic boundary condition,
# we only use a subset of the longitudes:
lon = ncread(fhret,"longitude")
indu = find( x -> x < 360., lon )
lon = lon[indu]
lat = ncread(fhret,"latitude")

println("Load harmonic constants and compute interpolators for all tides.")
k = 0
for cid in cidvec
    println("Loading")
    @time begin
        re  = ncread(fhret,cid * "re",start=[1,1],count=[length(indu),-1])
        im  = ncread(fhret,cid * "im",start=[1,1],count=[length(indu),-1])
    end

    if (FIG)
        cmap = Colors.colormap("RdBu")
        cmap = reverse(cmap)
        Winston.colormap(cmap)
        figure(height=400,width=800)
        p1=imagesc((lon[1],lon[end]),(lat[end],lat[1]),re',(-0.02,0.02))
        display(p1)
    end
    
    R = 6.3781e6 ; # meters
    n,m = size(re)
    
    dlon = lon[2] - lon[1]
    dlat = lat[2] - lat[1]
    
    knots = (lon,lat)
    zz = re + complex(0.0,1.0).*im
    itp=interpolate(zz,(
                        BSpline(Quadratic(Periodic(OnGrid()))),
                        BSpline(Quadratic(Reflect(OnGrid())))
                        ))
    
    zz_x = Array{Complex,2}(undef,size(zz))
    zz_y = Array{Complex,2}(undef,size(zz))
    @time begin
        for j=1:m
            println("j = ",j)
            for i=1:n
                dx = dlon*R*pi/180.0*cos( lat[j]*pi/180. )
                dy = dlat*R*pi/180.0
                ri = convert(Float64,i)
#            ri = collect(1.0:1.0*n)
                rj = convert(Float64,j)
                g = Interpolations.gradient(itp,ri,rj)
                zz_x[i,j] = g[1]./dx
                zz_y[i,j] = g[2]./dy
            end
        end
    end
    
    # This seems to work fine!
    # The components of the complex gradient are in zz_x and zz_y
    sm1 = complex(0.0,1.0)
    dtr = pi/180.
    rot = 2.0*pi/86164.1
    cori = zeros(size(lat))
    # for the grav0 solution:
    grav = ones(size(lat)).*9.80665
    for j=collect(1:m)
        cori[j] = 2.0*rot.*sin( lat[j]*dtr )
        # This is the International Gravity Formula, according to the Wikipedia:
        grav[j] = 9.780327*(1. + 0.0053024*sin(lat[j]*dtr).^2 - 0.0000058*sin(2.0*lat[j]*dtr).^2)
    end
    # for the grav0 solution:
    #grav = ones(size(lat)).*9.80665
    
    omega = 2.0*pi/(HAmod.period(cid)*8.64e4)
    
    println("Computing the velocity")
    xu = complex(zeros(n+1,m))
    xv = complex(zeros(n+1,m))
    for j=collect(1:m)
        J = complex(zeros(2,2))
        J[1,:] = [-sm1*omega+gamma,   -cori[j] ]
        J[2,:] = [   cori[j], -sm1*omega+gamma ]
        Jinv = inv(J)
        for i=collect(1:n)
            uv = -grav[j]*Jinv*[zz_x[i,j], zz_y[i,j]]
            xu[i,j] = uv[1]
            xv[i,j] = uv[2]
        end
        # Explicitly wrap to make periodic for linear interpolation:
        xu[n+1,j] = xu[1,j]
        xv[n+1,j] = xv[1,j]
    end
    
    spd = sqrt.( 0.5*(abs.(xu).^2 + abs.(xv).^2) )
    if (FIG)
        figure(height=400,width=800)
        p1=imagesc((lon[1],lon[end]),(lat[end],lat[1]),spd',(0.,0.1))
        display(p1)
    end
    
    println("Building velocity interpolators")
    @time begin
        # Now make interpolators for complex u and v:
        # Extend longitude by one unit for periodicity:
        lonex = [lon; lon[end]+dlon]
        knots=(lonex,lat)
        itpu=interpolate(knots,xu,(Gridded(Linear()),Gridded(Linear())))
        itpv=interpolate(knots,xv,(Gridded(Linear()),Gridded(Linear())))
    end

    global k = k + 1
    M[k] = Dict("cid"  => cid,
                "itpu" => itpu,
                "itpv" => itpv)
end

# load depth and distance from land
# build interpolators for depth and dist
lon_lim=[0. 361.]
lat_lim=[-66. 66.]
fmap = "/home/ezaron/NASA-Tides/GEBCO/GEBCO_gridone_distance_to_coastline_v2.0m.nc"
glon = nc_varget(fmap,"longitude")
glat = nc_varget(fmap,"latitude")
indi = find( x -> (x >= lon_lim[1]) & (x <= lon_lim[2]),glon)
indj = find( x -> (x >= lat_lim[1]) & (x <= lat_lim[2]),glat)
hght = nc_varget(fmap,"height",start=[indi[1],indj[1]],count=[length(indi),length(indj)])
glon = glon[indi]
glat = glat[indj]

# I should accumulate rms stats in a separate program, in case I want to 
# use specific frequencies or include a barotropic current:

# Set up arrays to accumulate the prior variance of velocity components
# and the posterior variance and the kount of drifter obs.
DLON=2.
DLAT=2.
# cell edges for collecting statistics:
slon=linspace( 38., 90. - DLON,round(Int,(90. - 38.)/DLON))
slat=linspace(-16., 24. - DLON,round(Int,(24. + 16.)/DLAT))
nn=length(slon)
mm=length(slat)
cnt = zeros(nn,mm)

ncid = length(M)
# Prior
ubarp=zeros(nn,mm)
vbarp=zeros(nn,mm)
uvarp=zeros(nn,mm)
vvarp=zeros(nn,mm)
# Tidal
ubart=zeros(nn,mm,ncid+1)
vbart=zeros(nn,mm,ncid+1)
uvart=zeros(nn,mm,ncid+1)
vvart=zeros(nn,mm,ncid+1)
# Final
ubarf=zeros(nn,mm,ncid+1)
vbarf=zeros(nn,mm,ncid+1)
uvarf=zeros(nn,mm,ncid+1)
vvarf=zeros(nn,mm,ncid+1)

fr = match(r"(.+?)\.nc",fu)
if ( fr == nothing )
    println("Error: drifters_predict.jl was not able to identify root of filename.")
    println("       Output will be written to the file tmp.nc.")
    fr = "tmp"
else
    fr = fr[1]
end
fn = fdir * fu

nci = ncinfo(fn)
xlon = ncread(fn,"LON")
xlat = ncread(fn,"LAT")
tau  = ncread(fn,"TIME")
u    = ncread(fn,"U")
v    = ncread(fn,"V")
driftid = ncread(fn,"ID")

# locate the ends of drifter trajectories:
# From the netcdf: "Interruptions in contiguous hourly estimates are indicated by \"Infinity\" values; 
# Individual trajectories are separated by \"NaN\" values"
indx = find(x -> x == Inf,xlon)
println("Identified ",length(indx)," uninterrupted hourly trajectories in fn = ",fn)
iend = indx .- 1
istart = [1 ; indx[1:end-1] .+ 1]
# Note that x == NaN always tests false! Must use isnan(x) instead.
indxx = find(isnan, xlon[istart])
println("Identified ",length(indxx)," (possibly-interrupted) drifter trajectories fn = ",fn)
istart[indxx] = istart[indxx] .+ 1

# convert time to jd:
tjd = HAmod.ymd2jd(1979.0,1.0,1.0) .+ tau.*24.

indx = find( x -> x < 0.,xlon )
xlon[indx] = xlon[indx] .+ 360.


############## SETUP THE OUTPUT FILE ###############

# We shall write a file which contains redundant time, lon, lat for matching with the
# original file, and which contains nan and inf markers in the same entries.
nd = length(tau)

# Create a netcdf file to hold the drifter data AND the tidal predictions:
fxout = fdir * fr * "_pred.nc"
# If the file already exists, then compute requested predictions and
# add these to the file. Otherwise make a fresh file.
fnsize = Base.Filesystem.filesize(fxout)
if (fnsize == 0)
    # Create a new file.
    # Everything is a function of the TIME coordinate dimension:
    nccreate(fxout,"tjd",
                      "TIME",tau,Dict("standard_name" => "time",
                                      "units" => "hours since 1979-01-01 00:00:00"),
             atts=Dict("standard_name" => "Julian date",
                       "units" => "real-valued Julian date")
             )
    nccreate(fxout,"LON",
             "TIME",
             atts=Dict("standard_name" => "longitude",
                       "units" => "deg E"))
    nccreate(fxout,"LAT",
             "TIME",
             atts=Dict("standard name" => "latitude",
                       "units" => "deg N"))
else
    # File already exists, append new predictions.
    ttmp = ncread(fxout,"TIME")
    if ( length(ttmp) != nd )
        println("Error: drifters_predict.jl: Output file ",fxout,
                " already exists; however, it appears to be the wrong size to append data.")
        "_break_"
    end
end
ncsync()

for cid=cidvec
    try
        nccreate(fxout,"u_hret_" * cid,
                 "TIME",
                 atts=Dict("standard_name" => "predicted " * cid * " u-vel", "units" => "m/s"))
        nccreate(fxout,"v_hret_" * cid,
                 "TIME",
                 atts=Dict("standard_name" => "predicted " * cid * " v-vel", "units" => "m/s"))
    catch
        println("Currents for cid = " * cid * " already exist in fxout = " * fxout)
    end
end

if (fnsize == 0)
    # Not sure what will happen if I attempt to write global attributes if they already exist in
    # the file:
    gatts = Dict("scripts"       => "namako:NASA-Tides/Julia/drifters_predict.jl",
                 "HRET_file"     => fhret,
                 "mask"          => "from the HRET file; see maskIZ.jl",
                 "drifter_file"  => fn,
                 "grav"          => "latitude-dependent International Gravity Formula",
                 "date"          => "2019-05-21",
                 "creator"       => "Ed Zaron, ezaron@pdx.edu",
                 "note"          => "Missing data, marked as nan or inf, should match exactly the drifter_file."
                 )
    ncputatt(fxout,"global",gatts)
    
    println("Writing original data to " * fu)
    ncwrite(xlon,fxout,"LON")
    ncwrite(xlat,fxout,"LAT")
    ncwrite(tau,fxout,"TIME")
    ncwrite(tjd,fxout,"tjd")
end

# initialize velocities with NaN's since some fall out of valid
# latitude range and so-forth.
nans = zeros(nd,).*NaN
for cid=cidvec
    ncwrite(nans,fxout,"u_hret_" * cid)
    ncwrite(nans,fxout,"v_hret_" * cid)
end

ncsync()

# See ../Julia_aji/pitvel.jl for example of tidal predictions.
k=300
#for k=collect(10000:11670)
for k=collect(1:length(istart))
    println("Working on set k=",k)
    ii=collect(istart[k]:iend[k])
    t1 = tjd[ii]
    lat1 = xlat[ii]
    lon1 = xlon[ii]

    indu = find( x ->
                (x >= lat[1]) &
                (x <= lat[end]),
                lat1 )

    if (length(indu) == 0)
        continue
    end

    ii = ii[indu]
    t1 = t1[indu]
    lat1 = lat1[indu]
    lon1 = lon1[indu]

    ic = 0
    for cid in cidvec
        # Nodal correction is turned on. Might want to figure out how to optimize this!
        # It might be possible to save time by building this matrix all at once
        # for all the cid's:
        iid,F = HAmod.FMAT(t1,[cid],1,0)

        # Harmonic constants:
        u0 = zeros(2,length(ii))
        v0 = zeros(2,length(ii))
        idc = iid[cid * "c"]
        ids = iid[cid * "s"]
        ic=ic+1
        itpu = M[ic]["itpu"]
        for i=collect(1:length(ii))
            u1 = itpu(lon1[i],lat1[i])
            u0[idc,i] = real(u1)
            u0[ids,i] = imag(u1)
        end
        
        itpv = M[ic]["itpv"]
        for i=collect(1:length(ii))
            v1 = itpv(lon1[i],lat1[i])
            v0[idc,i] = real(v1)
            v0[ids,i] = imag(v1)
        end

        F = transpose(F)
        u1p = sum(F.*u0,dims=1)
        v1p = sum(F.*v0,dims=1)

        ncwrite(u1p,fxout,"u_hret_" * cid,start=[ii[1]],count=[length(ii)])
        ncwrite(v1p,fxout,"v_hret_" * cid,start=[ii[1]],count=[length(ii)])
    end
    ncsync()
end

ncclose()