!*nictec -- Evaluate NIC at given location and time
!+
! These FORTRAN90 routines can be used to calculate the vertical total
! electron content based on the NOAA Ionosphere Climatology 2008 (NIC08).
! Later versions of the NIC climatology may be supported.
!
! The routines and a brief description:
! nictec  : Evaluate NIC at a given location and time
! nicinit : Initialize by reading the climatology and forcing files
! nicfree : Free up allocated memory
!
! To get some of the internal variables out of these routines, the user
! can insert the module nicclim into their code.
!
!-----------------------------------------------------------------------
! FUNCTION nictec (time, lat, lon)
! REAL(eightbytereal), INTENT(in) :: time, lat, lon
! REAL(eightbytereal) :: nictec
!
! This function determines the vertical total electron content at a given
! time (time) and location (lat, lon) based on the NOAA Ionosphere Climatology
! (NIC08) or later.
!
! The nictec routine needs to be initilized using a call to nicinit. This
! routine will load the NIC climatology file and a regularly updated file
! containing the forcing parameters for this climatology: 2-hourly values
! of the global mean TEC.
!
! The TEC is integrated to the height mentioned in the climatology files.
! Generally, two climatology files will be available (e.g. nic08_heo.nc and
! nic08_leo.nc) for high and low earth orbiter, respectively. See nicinit.
!
! The TEC value can be converted to ionospheric path delay by using the
! relationship:
!   PD = C * TEC / f**2
! where PD is path delay in metres, C is a constant of 40250d13 and f is the
! signal frequency in Hz.
!
! Input arguments are:
!  time    : time in UTC seconds since 1.0 Jan 1985
!  lat     : latitude in degrees
!  lon     : longitude in degrees (range from -540 to 540 is allowed)
!
! Returned value:
!  nictec  : total electron content (TEC units = 10^16 electrons/m^2)
!
!-----------------------------------------------------------------------
! SUBROUTINE nicinit (nicfile, gtecfile)
! CHARACTER(*), INTENT(in) :: nicfile, gtecfile
!
! This routine needs to be called prior to nictec. It reads the climatology
! file and the forcing parameter file into memory and performs some other
! initializations.
!
! The climatology files and forcing parameter file can be found at
! ftp://falcon.grdl.noaa.gov/pub/remko/nic08
!
! The nictec routine will halt the processing when the requested time is out of
! reach of the forcing parameter file (gtecfile). Hence that files should be
! updated reguraly.
! To allow some extrapolation, 7 days of GTEC values are automatically added
! to the end of the gtecfile.
!
! When nictec is no longer needed in your program, the memory allocated
! by nicinit can be freed using the nicfree routine.
!
! Input arguments are:
!  nicfile : path name of the NIC climatology file
!  gtecfile: path name of the GTEC forcing parameter time series
!
!-----------------------------------------------------------------------
! SUBROUTINE nicfree ()
!
! Free up the memory previously allocated by nicinit.
! The nicfree routine has no arguments.
!
!-----------------------------------------------------------------------
! MODULE nicclim
!
! To use this module as the line "use nicclim" to your code. However,
! this is absolutely not essential, unless you want to use some of the
! NIC's internal variables.
!
! nic_year        : First and last year of JPL GIM data used in NIC
! nic_gtec_tab    : Table of global mean TEC forcing
! nic_gtec_scale  : Scale factor to TECU
! nic_nt          : Number of entries in the GTEC table
! nic_trange      : Time range in the GTEC table (hours since 2000.0)
! nic_tec         : Array of the NIC climatology data
! nic_tec_scale   : Scale factor to TECU
! nic_level       : Lower and upper level of GTEC
! nic_gtec        : GTEC value at current time (after nictec call)
! nic_dtec        : Derivative of TEC to GTEC (after nictec call)
! 
!-----------------------------------------------------------------------
! $Log: nictec.f90,v $
! Revision 1.3  2008/01/25 21:27:16  rads
! - New version for NIC08. Now includes nicinit and nicfree
!
! Revision 1.1  2007/04/20 14:52:57  remko
! Initial revision
!
! Copyright (c) Remko Scharroo, Altimetrics LLC
!-----------------------------------------------------------------------

module nicclim
use typesizes
real(eightbytereal), save :: nic_tec_scale,nic_gtec_scale,nic_level(2),nic_trange(2),nic_gtec,nic_dtec
integer(fourbyteint), save :: nic_nt=0,nic_year(2)
integer(twobyteint), save, allocatable :: nic_gtec_tab(:),nic_tec(:,:,:,:,:)
end module

!-----------------------------------------------------------------------
! nicinit: Load NIC climatology and GTEC forcing file
!-----------------------------------------------------------------------

subroutine nicinit (nicfile, gtecfile)
use nicclim
use netcdf
character(*), intent(in) :: nicfile,gtecfile
integer(fourbyteint) :: ncid,varid
character(32) :: text
real(eightbytereal) :: gtec

if (nic_nt > 0) call nicfin("nicinit: NIC already initialized. Use nicfree first.")

! Read NIC climatology

call nfs(nf90_open(nicfile,nf90_nowrite,ncid))
call nfs(nf90_get_att(ncid,nf90_global,"description",text))
if (text(:5)=="NIC07") call nicfin("nicinit: This version does not work with NIC07")
call nfs(nf90_get_att(ncid,nf90_global,"year_range",nic_year))
call nfs(nf90_inq_varid(ncid,"gtec",varid))
call nfs(nf90_get_var(ncid,varid,nic_level))
call nfs(nf90_inq_varid(ncid,"tec",varid))
call nfs(nf90_get_att(ncid,varid,"scale_factor",nic_tec_scale))
allocate (nic_tec(73,73,12,12,2))
call nfs(nf90_get_var(ncid,varid,nic_tec))
call nfs(nf90_close(ncid))

! Make nic_level(2) the GTEC range (to save some computation later on)

nic_level(2) = nic_level(2) - nic_level(1)

! If gtecfile starts with "gtec=", simply set a constant GTEC

if (gtecfile(:5) == "gtec=") then
    nic_nt=2
    nic_gtec_scale=1d-2
    nic_trange=(/-1d20,2d20/)
    allocate (nic_gtec_tab(nic_nt))
    read (gtecfile(6:),*) gtec
    nic_gtec_tab=nint(gtec/nic_gtec_scale,twobyteint)
    return
endif

! Read GTEC data

call nfs(nf90_open(gtecfile,nf90_nowrite,ncid))
call nfs(nf90_inquire_dimension(ncid,1,len=nic_nt))
allocate (nic_gtec_tab(nic_nt))
call nfs(nf90_inq_varid(ncid,"time",varid))
call nfs(nf90_get_att(ncid,varid,"actual_range",nic_trange))
call nfs(nf90_inq_varid(ncid,"gtec",varid))
call nfs(nf90_get_att(ncid,varid,"scale_factor",nic_gtec_scale))
call nfs(nf90_get_var(ncid,varid,nic_gtec_tab))
call nfs(nf90_close(ncid))

! Make nic_trange(2) the time step (to save some computation later on)

nic_trange(2) = (nic_trange(2) - nic_trange(1)) / (nic_nt - 1)

contains

subroutine nfs(ios)
! Error exit for NetCDF calls
integer, intent(in) :: ios
if (ios==nf90_noerr) return
call nicfin("nicinit: "//nf90_strerror(ios))
end subroutine nfs

end subroutine nicinit

!-----------------------------------------------------------------------
! nicfin: General exit for all NIC routines
!-----------------------------------------------------------------------

subroutine nicfin(string)
character(*), intent(in) :: string
character(80) :: prognm
call getarg(0,prognm)
write (0,'(a,": ",a)') trim(prognm),trim(string)
stop
end subroutine nicfin

!-----------------------------------------------------------------------
! nicfree: Free up memory claimed by nicinit
!-----------------------------------------------------------------------

subroutine nicfree ()
use nicclim
if (nic_nt > 0) deallocate (nic_gtec_tab,nic_tec)
nic_nt = 0
end subroutine nicfree

!-----------------------------------------------------------------------
! nictec: Interpolate NIC at given time, latitude and longitude
!-----------------------------------------------------------------------

function nictec (time, lat, lon)
use netcdf
use nicclim
real(eightbytereal) :: nictec
real(eightbytereal), intent(in) :: time,lat,lon
integer(fourbyteint) :: kt,h0,h1,m0,m1,kx,ky
real(eightbytereal) :: t,f,xtime,h,m,x,y,upy,tec_sqr(2,2)
real(eightbytereal), parameter :: rad = atan(1d0)/45d0, lonp=288.63d0

if (nic_nt <= 0) call nicfin("nictec: NIC not initialized. Use nicinit first.")

! Convert time (UTC85) into hours since 2000.0

xtime = (time - 473299200d0) / 3600d0

! Determine location of matching GTEC record

t = (xtime - nic_trange(1)) / nic_trange(2) + 1d0
kt = floor(t)
t = t - kt
if (kt < 1 .or. kt >= nic_nt) call nicfin ("nictec: Time outside GTEC file")

! Interpolate forcing parameter (GTEC) and determine its fraction.

nic_gtec = (nic_gtec_tab(kt) * (1d0 - t) + nic_gtec_tab(kt+1) * t) * nic_gtec_scale
f = (nic_gtec - nic_level(1)) / nic_level(2)

! Determine month and month fraction. m0 and m1 are neighbouring months.

m = modulo (xtime / 730.5d0 + 0.5d0, 12d0)
m0 = floor(m)
m = m - m0
m1 = m0 + 1
if (m0 == 0) m0 = 12

! Determine bihour and bihour fraction. h0 and h1 are neighbouring bihours.

h = modulo (xtime / 2d0, 12d0) + 1d0
h0 = floor(h)
h = h - h0
h1 = h0 + 1
if (h1 == 13) h1 = 1

! Determine X coordinate, shifted East for earlier grid, and X fraction

x = modulo(lon + h * 30d0 + 180d0, 360d0) / 5d0 + 1d0
kx = floor(x)
x = x - kx

! Determine Y coordinate and Y fraction of grid point in the earlier grid

if (h >= 1d-7) then
    upy = 30d0 * 0.2d0 * sin((lon-lonp)*rad) * cos(lat*rad)
else	! Save a bit of CPU when h is tiny
    upy = 0d0
endif
y = (lat + h * upy + 90d0) / 2.5d0 + 1d0
ky = max(1,min(floor(y),72))
y = y - ky

! write (*,*) 't=',kt,kt+1,t
! write (*,*) 'x=',kx,kx+1,x
! write (*,*) 'y=',ky,ky+1,y
! write (*,*) 'm=',m0,m1,m
! write (*,*) 'h=',h0,h1,h
! write (*,*) 'f=',1,2,f

! Store earlier grid data while interpolating in latitude, month and
! bihour, so that tec_sqr will be a 2x2 matrix of which corners are the
! neighbouring longitude and forcing level.

tec_sqr = (1d0-h) * &
    ((1d0-y) * ((1d0-m) * nic_tec(kx:kx+1,ky  ,h0,m0,:) + m * nic_tec(kx:kx+1,ky  ,h0,m1,:)) &
        + y  * ((1d0-m) * nic_tec(kx:kx+1,ky+1,h0,m0,:) + m * nic_tec(kx:kx+1,ky+1,h0,m1,:)))

! Next part is for interpolating the later grid. This can be skipped when
! h is tiny

if (h >= 1d-7) then

! Determine X coordinate, shifted West for later grid (fraction stays the same).

kx = kx - 6
if (kx < 1) kx = kx + 72

! Determine Y coordinate and Y fraction of grid point in the later grid

y = (lat - (1d0-h) * upy + 90d0) / 2.5d0 + 1d0
ky = max(1,min(floor(y),72))
y = y - ky

! Add later grid data while interpolating in month and weighting by bihour

tec_sqr = tec_sqr + h * &
    ((1d0-y) * ((1d0-m) * nic_tec(kx:kx+1,ky  ,h1,m0,:) + m * nic_tec(kx:kx+1,ky  ,h1,m1,:)) &
        + y  * ((1d0-m) * nic_tec(kx:kx+1,ky+1,h1,m0,:) + m * nic_tec(kx:kx+1,ky+1,h1,m1,:)))

endif

! Interpolate in the remaining dimensions (longitude and forcing)

tec_sqr(1,:) = (1d0-x) * tec_sqr(1,:) + x * tec_sqr(2,:)
nictec = ((1d0-f) * tec_sqr(1,1) + f * tec_sqr(1,2)) * nic_tec_scale

! Compute the derivative of TEC to GTEC

nic_dtec = (tec_sqr(1,2) - tec_sqr(1,1)) * nic_tec_scale / nic_level(2)

end function nictec