MODULE io_netcdf
!-----------------------------------------------------------------------
! NAME
!     io_netcdf
!
! DESCRIPTION
!     Provides input/output procedures and parameters for the PEM.
!
! AUTHORS & DATE
!     JB Clement, 12/2025
!
! NOTES
!
!-----------------------------------------------------------------------

! DEPEDENCIES
! -----------
use numerics, only: dp, di, k4
use netcdf,   only: nf90_double, nf90_noerr, nf90_strerror, nf90_write, nf90_nowrite,                &
                    nf90_open, nf90_close, nf90_redef, nf90_enddef, nf90_inquire, nf90_max_var_dims, &
                    nf90_inq_dimid, nf90_inquire_dimension, nf90_inq_varid, nf90_inquire_variable,   &
                    nf90_def_var, nf90_get_var, nf90_put_var, nf90_get_att, nf90_put_att
use stoppage, only: stop_clean

! DECLARATION
! -----------
implicit none

! PARAMETERS
! ----------
character(8),  parameter :: start_name     = 'start.nc'
character(11), parameter :: start1D_name   = 'start1D.txt'
character(10), parameter :: startfi_name   = 'startfi.nc'
character(11), parameter :: startpem_name  = 'startpem.nc'
character(19), parameter :: xios_day_name1 = 'Xoutdaily4pem_Y1.nc'  ! XIOS daily output file, second to last PCM year
character(19), parameter :: xios_day_name2 = 'Xoutdaily4pem_Y2.nc'  ! XIOS daily output file, last PCM year
character(20), parameter :: xios_yr_name1  = 'Xoutyearly4pem_Y1.nc' ! XIOS yearly output file, second to last PCM year
character(20), parameter :: xios_yr_name2  = 'Xoutyearly4pem_Y2.nc' ! XIOS yearly output file, last PCM year
character(11), parameter :: diagevol_name  = 'diagevol.nc'

! VARIABLES
! ---------
logical(k4), private :: open_ncfile = .false.   ! Flag to true if a NetCDF file is already open
logical(k4), private :: open_diagevol = .false. ! Flag to true if "diagevol.nc" is already open
integer(di), private :: ncid_file               ! File ID
integer(di), private :: varid                   ! Variable ID
integer(di)          :: ncid_diagevol           ! File ID specific to "diagevol.nc"

! INTERFACES
! ----------
interface get_var_nc
    module procedure get_var0d_nc, get_var1d_nc, get_var2d_nc, get_var3d_nc, get_var4d_nc
end interface get_var_nc

interface put_var_nc
    module procedure put_var0d_nc, put_var1d_nc, put_var2d_nc, put_var3d_nc, put_var4d_nc
end interface put_var_nc

interface check_valid_var_nc
    module procedure check_valid_var0d_nc, check_valid_var1d_nc, check_valid_var2d_nc, check_valid_var3d_nc, check_valid_var4d_nc
end interface check_valid_var_nc

contains
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

!=======================================================================
SUBROUTINE check_nc(ierr,msg,found)
!-----------------------------------------------------------------------
! NAME
!     check_nc
!
! DESCRIPTION
!     NetCDF error handler.
!
! AUTHORS & DATE
!     JB Clement, 12/2025
!
! NOTES
!
!-----------------------------------------------------------------------

! DEPENDENCIES
! ------------
use display, only: print_err

! ARGUMENTS
! ---------
integer(di),  intent(in)            :: ierr
character(*), intent(in)            :: msg
logical(k4),  intent(out), optional :: found

! LOCAL VARIABLES
! ---------------
logical(k4) :: tmp_found

! CODE
! ----
if (ierr == nf90_noerr) then
    tmp_found = .true.
else
    tmp_found = .false.
end if

if (present(found)) then
    found = tmp_found
else
    if (.not. tmp_found) then
        call print_err(trim(nf90_strerror(ierr)))
        call stop_clean(__FILE__,__LINE__,'NetCDF error when '//trim(msg),ierr)
    end if
end if

END SUBROUTINE check_nc
!=======================================================================

!=======================================================================
SUBROUTINE open_nc(filename,mode,itime)
!-----------------------------------------------------------------------
! NAME
!     open_nc
!
! DESCRIPTION
!     Open a netCDF file for reading.

! AUTHORS & DATE
!     JB Clement, 12/2025

! NOTES
!
!-----------------------------------------------------------------------

! DECLARATION
! -----------
implicit none

! ARGUMENTS
! ---------
character(*), intent(in)            :: filename, mode
integer(di),  intent(out), optional :: itime

! CODE
! ----
! Diagevol logic
if (adjustl(trim(filename)) == diagevol_name) then
    if (adjustl(trim(mode)) == 'read') then
        call stop_clean(__FILE__,__LINE__,'opening mode "read" cannot be used with '//trim(filename)//'"!',1)
    else if (adjustl(trim(mode)) == 'write') then
        if (.not. open_diagevol) then
            call check_nc(nf90_open(trim(filename),nf90_write,ncid_diagevol),'opening file '//trim(filename)//' to write')
            if (present(itime)) call get_next_itime_nc('Time',itime) ! Next time index
            open_diagevol = .true.
        end if
        return
    else
        call stop_clean(__FILE__,__LINE__,'opening mode "'//mode//'" unknown!',1)
    end if
end if

! Standard logic
if (open_ncfile) then ! A file is already opened
    call stop_clean(__FILE__,__LINE__,'a NetCDF file is already opened!',1)
else
    if (adjustl(trim(mode)) == 'read') then
        call check_nc(nf90_open(trim(filename),nf90_nowrite,ncid_file),'opening file '//trim(filename)//' to read')
    else if (adjustl(trim(mode)) == 'write') then
        call check_nc(nf90_open(trim(filename),nf90_write,ncid_file),'opening file '//trim(filename)//' to write')
        if (present(itime)) call get_next_itime_nc('Time',itime) ! Next time index
    else
        call stop_clean(__FILE__,__LINE__,'opening mode "'//mode//'" unknown!',1)
    end if
    open_ncfile = .true.
end if

END SUBROUTINE open_nc
!=======================================================================

!=======================================================================
SUBROUTINE close_nc(filename)
!-----------------------------------------------------------------------
! NAME
!     close_nc
!
! DESCRIPTION
!     Open a netCDF file.

! AUTHORS & DATE
!     JB Clement, 12/2025

! NOTES
!
!-----------------------------------------------------------------------

! DECLARATION
! -----------
implicit none

! ARGUMENTS
! ---------
character(*), intent(in) :: filename

! CODE
! ----
if (adjustl(trim(filename)) == diagevol_name) then ! Diagevol logic
    call check_nc(nf90_close(ncid_diagevol),'closing file '//trim(filename))
    open_diagevol = .false.
else ! Standard logic
    call check_nc(nf90_close(ncid_file),'closing file '//trim(filename))
    open_ncfile = .false.
end if

END SUBROUTINE close_nc
!=======================================================================

!=======================================================================
SUBROUTINE get_dim_nc(dim_name,d,found)
!-----------------------------------------------------------------------
! NAME
!     get_dim_nc
!
! DESCRIPTION
!     Read a 0D variable from open NetCDF file.
!
! AUTHORS & DATE
!     JB Clement, 12/2025
!
! NOTES
!
!-----------------------------------------------------------------------

! DECLARATION
! -----------
implicit none

! ARGUMENTS
! ---------
character(*), intent(in)            :: dim_name ! Variable name
integer(di),  intent(out)           :: d        ! Output dimension
logical(k4),  intent(out), optional :: found

! LOCAL VARIABLES
! ---------------
integer(di) :: dimid ! Dimension ID

! CODE
! ----
if (present(found)) then
    call check_nc(nf90_inq_dimid(ncid_file,dim_name,dimid),'inquiring '//dim_name,found)
    call check_nc(nf90_inquire_dimension(ncid_file,dimid,len = d),'getting '//dim_name,found)
else
    call check_nc(nf90_inq_dimid(ncid_file,dim_name,dimid),'inquiring '//dim_name)
    call check_nc(nf90_inquire_dimension(ncid_file,dimid,len = d),'getting '//dim_name)
end if

END SUBROUTINE get_dim_nc
!=======================================================================

!=======================================================================
SUBROUTINE def_var_nc(var_name,title,units,dimids)
!-----------------------------------------------------------------------
! NAME
!     def_var_nc
!
! DESCRIPTION
!     Define a variable into open NetCDF file.
!
! AUTHORS & DATE
!     JB Clement, 01/2026
!
! NOTES
!
!-----------------------------------------------------------------------

! DECLARATION
! -----------
implicit none

! ARGUMENTS
! ---------
character(*),              intent(in) :: var_name, title, units ! Variable|title|units name
integer(di), dimension(:), intent(in) :: dimids                 ! Dimensions IDs

! LOCAL VARIABLES
! ---------------
logical(k4) :: found
integer(di) :: ncid

! CODE
! ----
! Diagevol logic priority over standard logic
if (open_diagevol) then
    ncid = ncid_diagevol
else
    ncid = ncid_file
end if

! Check if the variable exists
call check_nc(nf90_inq_varid(ncid,var_name,varid),'inquiring '//var_name,found)
if (found) return ! Variable is already defined

! Enter define mode
call check_nc(nf90_redef(ncid),'entering define mode')

! Define variable
call check_nc(nf90_def_var(ncid,var_name,NF90_DOUBLE,dimids,varid),'defining variable '//var_name)
call check_nc(nf90_put_att(ncid,varid,'title',title),'putting title attribute for '//var_name)
call check_nc(nf90_put_att(ncid,varid,'units',units),'putting units attribute for '//var_name)

! Leave define mode
call check_nc(nf90_enddef(ncid),'leaving define mode')

END SUBROUTINE def_var_nc
!=======================================================================

!=======================================================================
SUBROUTINE get_var0d_nc(var_name,var,found)
!-----------------------------------------------------------------------
! NAME
!     get_var0d_nc
!
! DESCRIPTION
!     Read a 0D variable from open NetCDF file.
!
! AUTHORS & DATE
!     JB Clement, 12/2025
!
! NOTES
!
!-----------------------------------------------------------------------

! DECLARATION
! -----------
implicit none

! ARGUMENTS
! ---------
character(*), intent(in)            :: var_name ! Variable name
real(dp),     intent(out)           :: var      ! Output variable
logical(k4),  intent(out), optional :: found

! CODE
! ----
if (present(found)) then
    call check_nc(nf90_inq_varid(ncid_file,var_name,varid),'inquiring '//var_name,found)
    call check_nc(nf90_get_var(ncid_file,varid,var),'getting '//var_name,found)
else
    call check_nc(nf90_inq_varid(ncid_file,var_name,varid),'inquiring '//var_name)
    call check_nc(nf90_get_var(ncid_file,varid,var),'getting '//var_name)
end if
call check_valid_var_nc(var_name,var)

END SUBROUTINE get_var0d_nc
!=======================================================================

!=======================================================================
SUBROUTINE get_var1d_nc(var_name,var,found)
!-----------------------------------------------------------------------
! NAME
!     get_var1d_nc
!
! DESCRIPTION
!     Read a 1D variable from open NetCDF file.
!
! AUTHORS & DATE
!     JB Clement, 12/2025
!
! NOTES
!
!-----------------------------------------------------------------------

! DECLARATION
! -----------
implicit none

! ARGUMENTS
! ---------
character(*),           intent(in)            :: var_name ! Variable name
real(dp), dimension(:), intent(out)           :: var      ! Output variable
logical(k4),            intent(out), optional :: found

! CODE
! ----
if (present(found)) then
    call check_nc(nf90_inq_varid(ncid_file,var_name,varid),'inquiring '//var_name,found)
    call check_nc(nf90_get_var(ncid_file,varid,var),'getting '//var_name,found)
else
    call check_nc(nf90_inq_varid(ncid_file,var_name,varid),'inquiring '//var_name)
    call check_nc(nf90_get_var(ncid_file,varid,var),'getting '//var_name)
end if
call check_valid_var_nc(var_name,var)

END SUBROUTINE get_var1d_nc
!=======================================================================

!=======================================================================
SUBROUTINE get_var2d_nc(var_name,var,found)
!-----------------------------------------------------------------------
! NAME
!     get_var2d_nc
!
! DESCRIPTION
!     Read a 2D variable from open NetCDF file.
!
! AUTHORS & DATE
!     JB Clement, 12/2025
!
! NOTES
!
!-----------------------------------------------------------------------

! DECLARATION
! -----------
implicit none

! ARGUMENTS
! ---------
character(*),             intent(in)            :: var_name ! Variable name
real(dp), dimension(:,:), intent(out)           :: var      ! Output variable
logical(k4),              intent(out), optional :: found

! CODE
! ----
if (present(found)) then
    call check_nc(nf90_inq_varid(ncid_file,var_name,varid),'inquiring '//var_name,found)
    call check_nc(nf90_get_var(ncid_file,varid,var),'getting '//var_name,found)
else
    call check_nc(nf90_inq_varid(ncid_file,var_name,varid),'inquiring '//var_name)
    call check_nc(nf90_get_var(ncid_file,varid,var),'getting '//var_name)
end if
call check_valid_var_nc(var_name,var)

END SUBROUTINE get_var2d_nc
!=======================================================================

!=======================================================================
SUBROUTINE get_var3d_nc(var_name,var,found)
!-----------------------------------------------------------------------
! NAME
!     get_var3d_nc
!
! DESCRIPTION
!     Read a 3D variable from open NetCDF file.
!
! AUTHORS & DATE
!     JB Clement, 12/2025
!
! NOTES
!
!-----------------------------------------------------------------------

! DECLARATION
! -----------
implicit none

! ARGUMENTS
! ---------
character(*),               intent(in)            :: var_name ! Variable name
real(dp), dimension(:,:,:), intent(out)           :: var      ! Output variable
logical(k4),                intent(out), optional :: found

! CODE
! ----
if (present(found)) then
    call check_nc(nf90_inq_varid(ncid_file,var_name,varid),'inquiring '//var_name,found)
    call check_nc(nf90_get_var(ncid_file,varid,var),'getting '//var_name,found)
else
    call check_nc(nf90_inq_varid(ncid_file,var_name,varid),'inquiring '//var_name)
    call check_nc(nf90_get_var(ncid_file,varid,var),'getting '//var_name)
end if
call check_valid_var_nc(var_name,var)

END SUBROUTINE get_var3d_nc
!=======================================================================

!=======================================================================
SUBROUTINE get_var4d_nc(var_name,var,found)
!-----------------------------------------------------------------------
! NAME
!     get_var4d_nc
!
! DESCRIPTION
!     Read a 4D variable from open NetCDF file.
!
! AUTHORS & DATE
!     JB Clement, 12/2025
!
! NOTES
!
!-----------------------------------------------------------------------

! DECLARATION
! -----------
implicit none

! ARGUMENTS
! ---------
character(*),                 intent(in)            :: var_name ! Variable name
real(dp), dimension(:,:,:,:), intent(out)           :: var      ! Output variable
logical(k4),                  intent(out), optional :: found

! CODE
! ----
if (present(found)) then
    call check_nc(nf90_inq_varid(ncid_file,var_name,varid),'inquiring '//var_name,found)
    call check_nc(nf90_get_var(ncid_file,varid,var),'getting '//var_name,found)
else
    call check_nc(nf90_inq_varid(ncid_file,var_name,varid),'inquiring '//var_name)
    call check_nc(nf90_get_var(ncid_file,varid,var),'getting '//var_name)
end if
call check_valid_var_nc(var_name,var)

END SUBROUTINE get_var4d_nc
!=======================================================================

!=======================================================================
SUBROUTINE put_var0d_nc(var_name,var,itime)
!-----------------------------------------------------------------------
! NAME
!     put_var0d_nc
!
! DESCRIPTION
!     Write a 0D variable into open NetCDF file.
!
! AUTHORS & DATE
!     JB Clement, 01/2026
!
! NOTES
!
!-----------------------------------------------------------------------

! DECLARATION
! -----------
implicit none

! ARGUMENTS
! ---------
character(*), intent(in)           :: var_name ! Variable name
real(dp),     intent(in)           :: var      ! Input variable
integer(di),  intent(in), optional :: itime    ! Current time index

! LOCAL VARIABLES
! ---------------
integer(di)                               :: ndims, unlimdimid, ncid
integer(di), dimension(nf90_max_var_dims) :: dimids
logical(k4)                               :: has_time

! CODE
! ----
! Diagevol logic priority over standard logic
if (open_diagevol) then
    ncid = ncid_diagevol
else
    ncid = ncid_file
end if

! Check if the variable holds a Time dimension (unlimited dim)
call check_nc(nf90_inq_varid(ncid,var_name,varid),'inquiring '//var_name)
call check_nc(nf90_inquire_variable(ncid,varid,ndims = ndims,dimids = dimids),'inquiring dims of '//var_name)
if (ndims > 0) then
    call check_nc(nf90_inquire(ncid,unlimitedDimId = unlimdimid),'inquiring unlimited dim')
    has_time = (dimids(ndims) == unlimdimid)
else
    has_time = .false.
end if

if (present(itime)) then ! Time-dependent write
    if (.not. has_time) call stop_clean(__FILE__,__LINE__,'itime present but variable has no Time dimension: '//var_name,1)
    ! For 0D variable with time, just write at the time index
    call check_nc(nf90_put_var(ncid,varid,var,start = (/itime/)),'putting '//var_name)
else ! Static write
    if (has_time) call stop_clean(__FILE__,__LINE__,'itime absent but variable is time-dependent: '//var_name,1)
    call check_nc(nf90_put_var(ncid,varid,var),'putting '//var_name)
end if

END SUBROUTINE put_var0d_nc
!=======================================================================

!=======================================================================
SUBROUTINE put_var1d_nc(var_name,var,itime)
!-----------------------------------------------------------------------
! NAME
!     put_var1d_nc
!
! DESCRIPTION
!     Write a 1D variable into open NetCDF file.
!
! AUTHORS & DATE
!     JB Clement, 01/2026
!
! NOTES
!
!-----------------------------------------------------------------------

! DECLARATION
! -----------
implicit none

! ARGUMENTS
! ---------
character(*),           intent(in)           :: var_name ! Variable name
real(dp), dimension(:), intent(in)           :: var      ! Input variable
integer(di),            intent(in), optional :: itime    ! Current time index

! LOCAL VARIABLES
! ---------------
integer(di)                               :: i, ndims, unlimdimid, ncid
integer(di), dimension(nf90_max_var_dims) :: dimids
integer(di), dimension(:), allocatable    :: strt, cnt
logical(k4)                               :: has_time

! CODE
! ----
! Diagevol logic priority over standard logic
if (open_diagevol) then
    ncid = ncid_diagevol
else
    ncid = ncid_file
end if

! Check if the variable holds a Time dimension (unlimited dim)
call check_nc(nf90_inq_varid(ncid,var_name,varid),'inquiring '//var_name)
call check_nc(nf90_inquire_variable(ncid,varid,ndims = ndims,dimids = dimids),'inquiring dims of '//var_name)
if (ndims > 0) then
    call check_nc(nf90_inquire(ncid,unlimitedDimId = unlimdimid),'inquiring unlimited dim')
    has_time = (dimids(ndims) == unlimdimid)
else
    has_time = .false.
end if

if (present(itime)) then ! Time-dependent write
    if (.not. has_time) call stop_clean(__FILE__,__LINE__,'itime present but variable has no Time dimension: '//var_name,1)
    allocate(strt(ndims),cnt(ndims))
    strt(:) = 1
    cnt(:) = 1
    strt(ndims) = itime
    do i = 1,ndims - 1
        cnt(i) = size(var,i)
    end do
    call check_nc(nf90_put_var(ncid,varid,var,start = strt,count = cnt),'putting '//var_name)
    deallocate(strt,cnt)
else ! Static write
    if (has_time) call stop_clean(__FILE__,__LINE__,'itime absent but variable is time-dependent: '//var_name,1)
    call check_nc(nf90_put_var(ncid,varid,var),'putting '//var_name)
end if

END SUBROUTINE put_var1d_nc
!=======================================================================

!=======================================================================
SUBROUTINE put_var2d_nc(var_name,var,itime)
!-----------------------------------------------------------------------
! NAME
!     put_var2d_nc
!
! DESCRIPTION
!     Write a 2D variable into open NetCDF file.
!
! AUTHORS & DATE
!     JB Clement, 01/2026
!
! NOTES
!
!-----------------------------------------------------------------------

! DECLARATION
! -----------
implicit none

! ARGUMENTS
! ---------
character(*),             intent(in)           :: var_name ! Variable name
real(dp), dimension(:,:), intent(in)           :: var      ! Input variable
integer(di),              intent(in), optional :: itime    ! Current time index

! LOCAL VARIABLES
! ---------------
integer(di)                               :: i, ndims, unlimdimid, ncid
integer(di), dimension(nf90_max_var_dims) :: dimids
integer(di), dimension(:), allocatable    :: strt, cnt
logical(k4)                               :: has_time

! CODE
! ----
! Diagevol logic priority over standard logic
if (open_diagevol) then
    ncid = ncid_diagevol
else
    ncid = ncid_file
end if

! Check if the variable holds a Time dimension (unlimited dim)
call check_nc(nf90_inq_varid(ncid,var_name,varid),'inquiring '//var_name)
call check_nc(nf90_inquire_variable(ncid,varid,ndims = ndims,dimids = dimids),'inquiring dims of '//var_name)
if (ndims > 0) then
    call check_nc(nf90_inquire(ncid,unlimitedDimId = unlimdimid),'inquiring unlimited dim')
    has_time = (dimids(ndims) == unlimdimid)
else
    has_time = .false.
end if

if (present(itime)) then ! Time-dependent write
    if (.not. has_time) call stop_clean(__FILE__,__LINE__,'itime present but variable has no Time dimension: '//var_name,1)
    allocate(strt(ndims),cnt(ndims))
    strt(:) = 1
    cnt(:) = 1
    strt(ndims) = itime
    do i = 1,ndims - 1
        cnt(i) = size(var,i)
    end do
    call check_nc(nf90_put_var(ncid,varid,var,start = strt,count = cnt),'putting '//var_name)
    deallocate(strt,cnt)
else ! Static write
    if (has_time) call stop_clean(__FILE__,__LINE__,'itime absent but variable is time-dependent: '//var_name,1)
    call check_nc(nf90_put_var(ncid,varid,var),'putting '//var_name)
end if

END SUBROUTINE put_var2d_nc
!=======================================================================

!=======================================================================
SUBROUTINE put_var3d_nc(var_name,var,itime)
!-----------------------------------------------------------------------
! NAME
!     put_var3d_nc
!
! DESCRIPTION
!     Write a 3D variable into open NetCDF file.
!
! AUTHORS & DATE
!     JB Clement, 01/2026
!
! NOTES
!
!-----------------------------------------------------------------------

! DECLARATION
! -----------
implicit none

! ARGUMENTS
! ---------
character(*),               intent(in)           :: var_name ! Variable name
real(dp), dimension(:,:,:), intent(in)           :: var      ! Input variable
integer(di),                intent(in), optional :: itime    ! Current time index

! LOCAL VARIABLES
! ---------------
integer(di)                               :: i, ndims, unlimdimid, ncid
integer(di), dimension(nf90_max_var_dims) :: dimids
integer(di), dimension(:), allocatable    :: strt, cnt
logical(k4)                               :: has_time

! CODE
! ----
! Diagevol logic priority over standard logic
if (open_diagevol) then
    ncid = ncid_diagevol
else
    ncid = ncid_file
end if

! Check if the variable holds a Time dimension (unlimited dim)
call check_nc(nf90_inq_varid(ncid,var_name,varid),'inquiring '//var_name)
call check_nc(nf90_inquire_variable(ncid,varid,ndims = ndims,dimids = dimids),'inquiring dims of '//var_name)
if (ndims > 0) then
    call check_nc(nf90_inquire(ncid,unlimitedDimId = unlimdimid),'inquiring unlimited dim')
    has_time = (dimids(ndims) == unlimdimid)
else
    has_time = .false.
end if

if (present(itime)) then ! Time-dependent write
    if (.not. has_time) call stop_clean(__FILE__,__LINE__,'itime present but variable has no Time dimension: '//var_name,1)
    allocate(strt(ndims),cnt(ndims))
    strt(:) = 1
    cnt(:) = 1
    strt(ndims) = itime
    do i = 1,ndims - 1
        cnt(i) = size(var,i)
    end do
    call check_nc(nf90_put_var(ncid,varid,var,start = strt,count = cnt),'putting '//var_name)
    deallocate(strt,cnt)
else ! Static write
    if (has_time) call stop_clean(__FILE__,__LINE__,'itime absent but variable is time-dependent: '//var_name,1)
    call check_nc(nf90_put_var(ncid,varid,var),'putting '//var_name)
end if

END SUBROUTINE put_var3d_nc
!=======================================================================

!=======================================================================
SUBROUTINE put_var4d_nc(var_name,var,itime)
!-----------------------------------------------------------------------
! NAME
!     put_var4d_nc
!
! DESCRIPTION
!     Write a 4D variable into open NetCDF file.
!
! AUTHORS & DATE
!     JB Clement, 01/2026
!
! NOTES
!
!-----------------------------------------------------------------------

! DECLARATION
! -----------
implicit none

! ARGUMENTS
! ---------
character(*),                 intent(in)           :: var_name ! Variable name
real(dp), dimension(:,:,:,:), intent(in)           :: var      ! Input variable
integer(di),                  intent(in), optional :: itime    ! Current time index

! LOCAL VARIABLES
! ---------------
integer(di)                               :: i, ndims, unlimdimid, ncid
integer(di), dimension(nf90_max_var_dims) :: dimids
integer(di), dimension(:), allocatable    :: strt, cnt
logical(k4)                               :: has_time

! CODE
! ----
! Diagevol logic priority over standard logic
if (open_diagevol) then
    ncid = ncid_diagevol
else
    ncid = ncid_file
end if

! Check if the variable holds a Time dimension (unlimited dim)
call check_nc(nf90_inq_varid(ncid,var_name,varid),'inquiring '//var_name)
call check_nc(nf90_inquire_variable(ncid,varid,ndims = ndims,dimids = dimids),'inquiring dims of '//var_name)
if (ndims > 0) then
    call check_nc(nf90_inquire(ncid,unlimitedDimId = unlimdimid),'inquiring unlimited dim')
    has_time = (dimids(ndims) == unlimdimid)
else
    has_time = .false.
end if

if (present(itime)) then ! Time-dependent write
    if (.not. has_time) call stop_clean(__FILE__,__LINE__,'itime present but variable has no Time dimension: '//var_name,1)
    allocate(strt(ndims),cnt(ndims))
    strt(:) = 1
    cnt(:) = 1
    strt(ndims) = itime
    do i = 1,ndims - 1
        cnt(i) = size(var,i)
    end do
    call check_nc(nf90_put_var(ncid,varid,var,start = strt,count = cnt),'putting '//var_name)
    deallocate(strt,cnt)
else ! Static write
    if (has_time) call stop_clean(__FILE__,__LINE__,'itime absent but variable is time-dependent: '//var_name,1)
    call check_nc(nf90_put_var(ncid,varid,var),'putting '//var_name)
end if

END SUBROUTINE put_var4d_nc
!=======================================================================

!=======================================================================
SUBROUTINE get_next_itime_nc(dim_name,itime)
!-----------------------------------------------------------------------
! NAME
!     get_next_itime_nc
!
! DESCRIPTION
!     Get the next time index in a NetCDF file to record variables.
!
! AUTHORS & DATE
!     JB Clement, 01/2026
!
! NOTES
!     In most cases, dim_name = 'Time'.
!-----------------------------------------------------------------------

! DECLARATION
! -----------
implicit none

! ARGUMENTS
! ---------
character(*), intent(in)  :: dim_name ! Dimension name
integer(di),  intent(out) :: itime    ! Time index

! LOCAL VARIABLES
! ---------------
integer(di) :: length
logical(k4) :: found

! CODE
! ----
! Get dimension length
call get_dim_nc(dim_name,length,found)
if (.not. found) call stop_clean(__FILE__,__LINE__,'dimension '//dim_name//' not found',1)

! Next time index
itime = length + 1

END SUBROUTINE get_next_itime_nc
!=======================================================================

!=======================================================================
SUBROUTINE check_valid_var0d_nc(var_name,var)
!-----------------------------------------------------------------------
! NAME
!     check_valid_var0d_nc
!
! DESCRIPTION
!     Check the validity of a 0D variable read from a NetCDF file.
!
! AUTHORS & DATE
!     JB Clement, 02/2026
!
! NOTES
!
!-----------------------------------------------------------------------

! DEPENDENCIES
! ------------
use numerics, only: largest_nb
use stoppage, only: stop_clean
use utility,  only: real2str

! DECLARATION
! -----------
implicit none

! ARGUMENTS
! ---------
character(*), intent(in) :: var_name ! Variable name
real(dp),     intent(in) :: var      ! Input variable

! LOCAL VARIABLES
! ---------------
logical(k4)            :: has_fill, has_range, has_valid_min, has_valid_max
real(dp)               :: fill_value, valid_min, valid_max
real(dp), dimension(2) :: valid_range
integer(di)            :: ncid

! CODE
! ----
! Diagevol logic priority over standard logic
if (open_diagevol) then
    ncid = ncid_diagevol
else
    ncid = ncid_file
end if

! NaN
if (var /= var) call stop_clean(__FILE__,__LINE__,'NaN detected in variable '//var_name//'!',1)

! Infinite
if (abs(var) > largest_nb) call stop_clean(__FILE__,__LINE__,'Infs detected in variable '//var_name//'!',1)

! Fill value
has_fill = .false.
call check_nc(nf90_get_att(ncid,varid,"_FillValue",fill_value),'getting fill value',has_fill)
if (.not. has_fill) call check_nc(nf90_get_att(ncid,varid,"missing_value",fill_value),'getting missing value',has_fill)
if (has_fill .and. var == fill_value) call stop_clean(__FILE__,__LINE__,'Fill values ('//real2str(fill_value)//') detected in '//var_name//'!',1)

! Valid range
has_range = .false.
call check_nc(nf90_get_att(ncid,varid,"valid_range",valid_range),'getting valid range',has_range)
if (has_range) then
    valid_min = valid_range(1)
    valid_max = valid_range(2)
else! valid_min / valid_max (fallback)
    has_valid_min = .false.
    has_valid_max = .false.
    call check_nc(nf90_get_att(ncid,varid,"valid_min",valid_min),'getting valid min',has_valid_min)
    call check_nc(nf90_get_att(ncid,varid,"valid_max",valid_max),'getting valid max',has_valid_max)
    if (has_valid_min .or. has_valid_max) then
        has_range = .true.
        if (.not. has_valid_min) valid_min = -largest_nb
        if (.not. has_valid_max) valid_max = largest_nb
    end if
end if
if (has_range .and. (var < valid_min .or. var > valid_max)) call stop_clean(__FILE__,__LINE__,'Values outside valid range ('//real2str(valid_min)//','//real2str(valid_max)//') detected in '//var_name//'!',1)

END SUBROUTINE check_valid_var0d_nc
!=======================================================================

!=======================================================================
SUBROUTINE check_valid_var1d_nc(var_name,var)
!-----------------------------------------------------------------------
! NAME
!     check_valid_var1d_nc
!
! DESCRIPTION
!     Check the validity of a 1D variable read from a NetCDF file.
!
! AUTHORS & DATE
!     JB Clement, 02/2026
!
! NOTES
!
!-----------------------------------------------------------------------

! DEPENDENCIES
! ------------
use numerics, only: largest_nb
use stoppage, only: stop_clean
use utility,  only: real2str

! DECLARATION
! -----------
implicit none

! ARGUMENTS
! ---------
character(*),           intent(in) :: var_name ! Variable name
real(dp), dimension(:), intent(in) :: var      ! Input variable

! LOCAL VARIABLES
! ---------------
logical(k4)            :: has_fill, has_range, has_valid_min, has_valid_max
real(dp)               :: fill_value, valid_min, valid_max
real(dp), dimension(2) :: valid_range
integer(di)            :: ncid

! CODE
! ----
! Diagevol logic priority over standard logic
if (open_diagevol) then
    ncid = ncid_diagevol
else
    ncid = ncid_file
end if

! NaN
if (any(var /= var)) call stop_clean(__FILE__,__LINE__,'NaN detected in variable '//var_name//'!',1)

! Infinite
if (any(abs(var) > largest_nb)) call stop_clean(__FILE__,__LINE__,'Infs detected in variable '//var_name//'!',1)

! Fill value
has_fill = .false.
call check_nc(nf90_get_att(ncid,varid,"_FillValue",fill_value),'getting fill value',has_fill)
if (.not. has_fill) call check_nc(nf90_get_att(ncid,varid,"missing_value",fill_value),'getting missing value',has_fill)
if (has_fill) then
    if (any(var == fill_value)) call stop_clean(__FILE__,__LINE__,'Fill values ('//real2str(fill_value)//') detected in '//var_name//'!',1)
end if

! Valid range
has_range = .false.
call check_nc(nf90_get_att(ncid,varid,"valid_range",valid_range),'getting valid range',has_range)
if (has_range) then
    valid_min = valid_range(1)
    valid_max = valid_range(2)
else! valid_min / valid_max (fallback)
    has_valid_min = .false.
    has_valid_max = .false.
    call check_nc(nf90_get_att(ncid,varid,"valid_min",valid_min),'getting valid min',has_valid_min)
    call check_nc(nf90_get_att(ncid,varid,"valid_max",valid_max),'getting valid max',has_valid_max)
    if (has_valid_min .or. has_valid_max) then
        has_range = .true.
        if (.not. has_valid_min) valid_min = -largest_nb
        if (.not. has_valid_max) valid_max = largest_nb
    end if
end if
if (has_range) then
    if (any(var < valid_min) .or. any(var > valid_max)) call stop_clean(__FILE__,__LINE__,'Values outside valid range ('//real2str(valid_min)//','//real2str(valid_max)//') detected in '//var_name//'!',1)
end if

END SUBROUTINE check_valid_var1d_nc
!=======================================================================

!=======================================================================
SUBROUTINE check_valid_var2d_nc(var_name,var)
!-----------------------------------------------------------------------
! NAME
!     check_valid_var2d_nc
!
! DESCRIPTION
!     Check the validity of a 1D variable read from a NetCDF file.
!
! AUTHORS & DATE
!     JB Clement, 02/2026
!
! NOTES
!
!-----------------------------------------------------------------------

! DEPENDENCIES
! ------------
use numerics, only: largest_nb
use stoppage, only: stop_clean
use utility,  only: real2str

! DECLARATION
! -----------
implicit none

! ARGUMENTS
! ---------
character(*),             intent(in) :: var_name ! Variable name
real(dp), dimension(:,:), intent(in) :: var      ! Input variable

! LOCAL VARIABLES
! ---------------
logical(k4)            :: has_fill, has_range, has_valid_min, has_valid_max
real(dp)               :: fill_value, valid_min, valid_max
real(dp), dimension(2) :: valid_range
integer(di)            :: ncid

! CODE
! ----
! Diagevol logic priority over standard logic
if (open_diagevol) then
    ncid = ncid_diagevol
else
    ncid = ncid_file
end if

! NaN
if (any(var /= var)) call stop_clean(__FILE__,__LINE__,'NaN detected in variable '//var_name//'!',1)

! Infinite
if (any(abs(var) > largest_nb)) call stop_clean(__FILE__,__LINE__,'Infs detected in variable '//var_name//'!',1)

! Fill value
has_fill = .false.
call check_nc(nf90_get_att(ncid,varid,"_FillValue",fill_value),'getting fill value',has_fill)
if (.not. has_fill) call check_nc(nf90_get_att(ncid,varid,"missing_value",fill_value),'getting missing value',has_fill)
if (has_fill .and. any(var == fill_value)) call stop_clean(__FILE__,__LINE__,'Fill values ('//real2str(fill_value)//') detected in '//var_name//'!',1)

! Valid range
has_range = .false.
call check_nc(nf90_get_att(ncid,varid,"valid_range",valid_range),'getting valid range',has_range)
if (has_range) then
    valid_min = valid_range(1)
    valid_max = valid_range(2)
else! valid_min / valid_max (fallback)
    has_valid_min = .false.
    has_valid_max = .false.
    call check_nc(nf90_get_att(ncid,varid,"valid_min",valid_min),'getting valid min',has_valid_min)
    call check_nc(nf90_get_att(ncid,varid,"valid_max",valid_max),'getting valid max',has_valid_max)
    if (has_valid_min .or. has_valid_max) then
        has_range = .true.
        if (.not. has_valid_min) valid_min = -largest_nb
        if (.not. has_valid_max) valid_max = largest_nb
    end if
end if
if (has_range .and. (any(var < valid_min) .or. any(var > valid_max))) call stop_clean(__FILE__,__LINE__,'Values outside valid range ('//real2str(valid_min)//','//real2str(valid_max)//') detected in '//var_name//'!',1)

END SUBROUTINE check_valid_var2d_nc
!=======================================================================

!=======================================================================
SUBROUTINE check_valid_var3d_nc(var_name,var)
!-----------------------------------------------------------------------
! NAME
!     check_valid_var3d_nc
!
! DESCRIPTION
!     Check the validity of a 1D variable read from a NetCDF file.
!
! AUTHORS & DATE
!     JB Clement, 02/2026
!
! NOTES
!
!-----------------------------------------------------------------------

! DEPENDENCIES
! ------------
use numerics, only: largest_nb
use stoppage, only: stop_clean
use utility,  only: real2str

! DECLARATION
! -----------
implicit none

! ARGUMENTS
! ---------
character(*),               intent(in) :: var_name ! Variable name
real(dp), dimension(:,:,:), intent(in) :: var      ! Input variable

! LOCAL VARIABLES
! ---------------
logical(k4)            :: has_fill, has_range, has_valid_min, has_valid_max
real(dp)               :: fill_value, valid_min, valid_max
real(dp), dimension(2) :: valid_range
integer(di)            :: ncid

! CODE
! ----
! Diagevol logic priority over standard logic
if (open_diagevol) then
    ncid = ncid_diagevol
else
    ncid = ncid_file
end if

! NaN
if (any(var /= var)) call stop_clean(__FILE__,__LINE__,'NaN detected in variable '//var_name//'!',1)

! Infinite
if (any(abs(var) > largest_nb)) call stop_clean(__FILE__,__LINE__,'Infs detected in variable '//var_name//'!',1)

! Fill value
has_fill = .false.
call check_nc(nf90_get_att(ncid,varid,"_FillValue",fill_value),'getting fill value',has_fill)
if (.not. has_fill) call check_nc(nf90_get_att(ncid,varid,"missing_value",fill_value),'getting missing value',has_fill)
if (has_fill .and. any(var == fill_value)) call stop_clean(__FILE__,__LINE__,'Fill values ('//real2str(fill_value)//') detected in '//var_name//'!',1)

! Valid range
has_range = .false.
call check_nc(nf90_get_att(ncid,varid,"valid_range",valid_range),'getting valid range',has_range)
if (has_range) then
    valid_min = valid_range(1)
    valid_max = valid_range(2)
else! valid_min / valid_max (fallback)
    has_valid_min = .false.
    has_valid_max = .false.
    call check_nc(nf90_get_att(ncid,varid,"valid_min",valid_min),'getting valid min',has_valid_min)
    call check_nc(nf90_get_att(ncid,varid,"valid_max",valid_max),'getting valid max',has_valid_max)
    if (has_valid_min .or. has_valid_max) then
        has_range = .true.
        if (.not. has_valid_min) valid_min = -largest_nb
        if (.not. has_valid_max) valid_max = largest_nb
    end if
end if
if (has_range .and. (any(var < valid_min) .or. any(var > valid_max))) call stop_clean(__FILE__,__LINE__,'Values outside valid range ('//real2str(valid_min)//','//real2str(valid_max)//') detected in '//var_name//'!',1)

END SUBROUTINE check_valid_var3d_nc
!=======================================================================

!=======================================================================
SUBROUTINE check_valid_var4d_nc(var_name,var)
!-----------------------------------------------------------------------
! NAME
!     check_valid_var4d_nc
!
! DESCRIPTION
!     Check the validity of a 1D variable read from a NetCDF file.
!
! AUTHORS & DATE
!     JB Clement, 02/2026
!
! NOTES
!
!-----------------------------------------------------------------------

! DEPENDENCIES
! ------------
use numerics, only: largest_nb
use stoppage, only: stop_clean
use utility,  only: real2str

! DECLARATION
! -----------
implicit none

! ARGUMENTS
! ---------
character(*),                 intent(in) :: var_name ! Variable name
real(dp), dimension(:,:,:,:), intent(in) :: var      ! Input variable

! LOCAL VARIABLES
! ---------------
logical(k4)            :: has_fill, has_range, has_valid_min, has_valid_max
real(dp)               :: fill_value, valid_min, valid_max
real(dp), dimension(2) :: valid_range
integer(di)            :: ncid

! CODE
! ----
! Diagevol logic priority over standard logic
if (open_diagevol) then
    ncid = ncid_diagevol
else
    ncid = ncid_file
end if

! NaN
if (any(var /= var)) call stop_clean(__FILE__,__LINE__,'NaN detected in variable '//var_name//'!',1)

! Infinite
if (any(abs(var) > largest_nb)) call stop_clean(__FILE__,__LINE__,'Infs detected in variable '//var_name//'!',1)

! Fill value
has_fill = .false.
call check_nc(nf90_get_att(ncid,varid,"_FillValue",fill_value),'getting fill value',has_fill)
if (.not. has_fill) call check_nc(nf90_get_att(ncid,varid,"missing_value",fill_value),'getting missing value',has_fill)
if (has_fill .and. any(var == fill_value)) call stop_clean(__FILE__,__LINE__,'Fill values ('//real2str(fill_value)//') detected in '//var_name//'!',1)

! Valid range
has_range = .false.
call check_nc(nf90_get_att(ncid,varid,"valid_range",valid_range),'getting valid range',has_range)
if (has_range) then
    valid_min = valid_range(1)
    valid_max = valid_range(2)
else! valid_min / valid_max (fallback)
    has_valid_min = .false.
    has_valid_max = .false.
    call check_nc(nf90_get_att(ncid,varid,"valid_min",valid_min),'getting valid min',has_valid_min)
    call check_nc(nf90_get_att(ncid,varid,"valid_max",valid_max),'getting valid max',has_valid_max)
    if (has_valid_min .or. has_valid_max) then
        has_range = .true.
        if (.not. has_valid_min) valid_min = -largest_nb
        if (.not. has_valid_max) valid_max = largest_nb
    end if
end if
if (has_range .and. (any(var < valid_min) .or. any(var > valid_max))) call stop_clean(__FILE__,__LINE__,'Values outside valid range ('//real2str(valid_min)//','//real2str(valid_max)//') detected in '//var_name//'!',1)

END SUBROUTINE check_valid_var4d_nc
!=======================================================================

END MODULE io_netcdf