program gev_fit
  implicit none
  
  integer, parameter :: n = 100   ! number of data points
  real, dimension(n) :: data      ! input data
  real :: loc, scale, shape       ! GEV parameters
  real :: loglik                  ! log-likelihood
  integer :: i                    ! loop index
  
  ! Read in the input data
  open(10, file='data.txt')
  do i = 1, n
    read(10,*) data(i)
  end do
  close(10)
  
  ! Estimate the GEV parameters using maximum likelihood estimation
  call gev_mle(data, n, loc, scale, shape, loglik)
  
  ! Output the results
  write(*,*) 'GEV parameters:'
  write(*,*) 'Location = ', loc
  write(*,*) 'Scale = ', scale
  write(*,*) 'Shape = ', shape
  write(*,*) 'Log-likelihood = ', loglik
  
contains

  ! Subroutine to estimate GEV parameters using maximum likelihood estimation
  subroutine gev_mle(data, n, loc, scale, shape, loglik)
    implicit none
    
    real, dimension(n) :: data  ! input data
    real :: loc, scale, shape   ! GEV parameters
    real :: loglik              ! log-likelihood
    integer :: n                ! number of data points
    real :: eps                 ! convergence tolerance
    integer :: i                ! loop index
    real :: u, x                ! transformed data
    
    ! Initialize parameters
    loc = mean(data)
    scale = std_dev(data)
    shape = 0.0
    
    ! Set convergence tolerance
    eps = 1e-6
    
    ! Iterate until convergence
    do while (abs(loglik - gev_loglik(data, n, loc, scale, shape)) > eps)
      loglik = gev_loglik(data, n, loc, scale, shape)
      shape = gev_shape(data, n, loc, scale)
      loc = gev_loc(data, n, scale, shape)
      scale = gev_scale(data, n, loc, shape)
    end do
    
    ! Calculate the final log-likelihood
    loglik = gev_loglik(data, n, loc, scale, shape)
    
  end subroutine gev_mle
  
  ! Function to calculate the log-likelihood of the GEV distribution
  real function gev_loglik(data, n, loc, scale, shape)
    implicit none
    
    real, dimension(n) :: data  ! input data
    real :: loc, scale, shape   ! GEV parameters
    integer :: n                ! number of data points
    real :: u, x                ! transformed data
    integer :: i                ! loop index
    
    ! Calculate the log-likelihood
    gev_loglik = 0.0
    do i = 1, n
      x = (data(i) - loc) / scale
      u = 1.0 + shape * x
      if (u > 0.0) then
        gev_loglik = gev_loglik - log(scale) - u**(-1.0/shape) - (1.0/shape+1.0) * log(u)
      else
        gev_loglik = -1.0e30
        return
      end if
    end do
    
  end function gev_loglik
  
! Function to estimate the shape parameter of the GEV distribution
real function gev_shape(data, n, loc, scale)
  implicit none
  
  real, dimension(n) :: data  ! input data
  real :: loc, scale          ! GEV parameters
  integer :: n                ! number of data points
  real :: u, x, y             ! transformed data
  integer :: i                ! loop index
  
  ! Calculate the sample mean and standard deviation of the transformed data
  x = (data - loc) / scale
  y = log(1.0 + x)
  u = sum(y) / n
  x = y - u
  u = sum(x**2) / (n-1)
  
  ! Calculate the shape parameter
  if (u <= 0.0) then
    gev_shape = 0.0
  else
    gev_shape = 0.5 / log(2.0) * log(u)
  end if
  
end function gev_shape

! Function to estimate the location parameter of the GEV distribution
real function gev_loc(data, n, scale, shape)
  implicit none
  
  real, dimension(n) :: data  ! input data
  real :: scale, shape        ! GEV parameters
  integer :: n                ! number of data points
  real :: u, x                ! transformed data
  integer :: i                ! loop index
  
  ! Calculate the location parameter
  if (shape == 0.0) then
    gev_loc = mean(data)
  else
    x = (data - mean(data)) / scale
    u = (1.0 - shape * x)**(-1.0/shape)
    gev_loc = mean(data) - scale * (u - 1.0) / shape
  end if
  
end function gev_loc

! Function to estimate the mean value of a dataset
real function mean(data, n)
  implicit none
  
  real, dimension(n) :: data  ! input data
  integer :: n                ! number of data points
  real :: mu                  ! sample mean
  integer :: i                ! loop index
  
  ! Calculate the sample mean
  mu = sum(data) / n
  
  ! Return the result
  mean = mu
  
end function mean

! Function to estimate the standard deviation of a dataset
real function std_dev(data, n)
  implicit none
  
  real, dimension(n) :: data  ! input data
  integer :: n                ! number of data points
  real :: mu, sigma           ! sample mean and standard deviation
  integer :: i                ! loop index
  
  ! Calculate the sample mean
  mu = sum(data) / n
  
  ! Calculate the sample standard deviation
  sigma = sqrt(sum((data - mu)**2) / (n-1))
  
  ! Return the result
  std_dev = sigma
  
end function std_dev

END