MODULE phys_constants
!-----------------------------------------------------------------------
! NAME
!     phys_constants
!
! DESCRIPTION
!     Physical constants used across PEM modules.
!
! AUTHORS & DATE
!     JB Clement, 12/2025
!
! NOTES
!     Constants are initialized from NetCDF 'startfi.nc' via read_constants.
!-----------------------------------------------------------------------

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

! MODULE PARAMETERS
! -----------------
real, parameter :: pi = 4.*atan(1.) ! PI = 3.14159...

real :: g     ! Gravity (m/s2)
real :: r     ! Reduced gas constant,r = 8.314511/(mugaz/1000.0)
real :: mugaz ! Molar mass of the atmosphere (g/mol)
real :: rad   ! Radius of the planet (m)
real :: cpp   ! Cp of the atmosphere
real :: rcp   ! r/cpp

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

!=======================================================================
SUBROUTINE read_constants(filename)
!-----------------------------------------------------------------------
! NAME
!     read_constants
!
! DESCRIPTION
!     Read physical constants from NetCDF file 'startfi.nc' in variable
!     'controle' and initialize module-level constants.
!
! AUTHORS & DATE
!     JB Clement, 12/2025
!
! NOTES
!     Reads controle(5,7,8,9) for rad, g, mugaz, rcp and computes r.
!-----------------------------------------------------------------------

! DEPENDENCIES
! ------------
use netcdf

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

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

! LOCAL VARIABLES
! ---------------
real, dimension(:), allocatable :: controle
integer :: ncid           ! File ID
integer :: varid_controle ! Variable ID for 'controle'
integer :: dimid_index    ! Dimension ID for 'index'
integer :: nindex         ! Size of dimension 'index'
integer :: ierr           ! Return codes

! CODE
! ----
! Open the NetCDF file
ierr = nf90_open(trim(filename),NF90_NOWRITE,ncid)
if (ierr /= nf90_noerr) then
    write(*,*) "Error opening file:", trim(nf90_strerror(ierr))
    error stop
endif

! Get the dimension size of 'index'
ierr = nf90_inq_dimid(ncid,"index",dimid_index)
if (ierr /= nf90_noerr) then
    write(*,*) "Error getting dimid 'index':", trim(nf90_strerror(ierr))
    error stop
endif
ierr = nf90_inquire_dimension(ncid,dimid_index,len = nindex)
if (ierr /= nf90_noerr) then
    write(*,*) "Error getting dimension length:", trim(nf90_strerror(ierr))
    error stop
endif

! Get the variable ID for 'controle'
allocate(controle(nindex))
ierr = nf90_inq_varid(ncid,"controle",varid_controle)
if (ierr /= nf90_noerr) then
    write(*,*) "Error getting variable ID 'controle':", trim(nf90_strerror(ierr))
    error stop
endif
ierr = nf90_get_var(ncid,varid_controle,controle)
if (ierr /= nf90_noerr) then
    write(*,*) "Error reading 'controle':", trim(nf90_strerror(ierr))
    error stop
endif

! Close the file
ierr = nf90_close(ncid)
if (ierr /= nf90_noerr) then
    write(*,*) "Error closing file:", trim(nf90_strerror(ierr))
    error stop
endif

! Initialize the constants with 'controle' data
rad   = controle(5)
g     = controle(7)
mugaz = controle(8)
rcp   = controle(9)
r     = 8.314511*1000./mugaz
cpp   = r/rcp
deallocate(controle)

END SUBROUTINE read_constants
!=======================================================================

END MODULE phys_constants