MODULE FRSEL_SETUP      ! specifies program parameters and sets-up haul information data type

  IMPLICIT NONE

  INTEGER, PARAMETER    :: nalpha = 4           ! number of parameters to be estimated
  INTEGER, PARAMETER	:: nnu = 2              ! number of selectivity parameters for each haul
  INTEGER, PARAMETER    :: wp = KIND(1.0d0)     ! kind-value for double precision

  TYPE HAUL_INFO                                ! derived type to contain haul information
    REAL (wp), DIMENSION(nnu,nalpha)  :: x      ! design matrix
    REAL (wp), DIMENSION(nnu)         :: nu, b  ! selectivity parameters and random effects
    REAL (wp), DIMENSION(nnu,nnu)     :: r, w   ! within-haul and (inverse of) total variance 
  END TYPE HAUL_INFO

END MODULE FRSEL_SETUP


INCLUDE 'FRSEL.F90'


PROGRAM EXAMPLE

  USE FRSEL_SETUP
  USE FRSEL

  ! example of the use of FRSEL1 taken from Fryer, 1991 (ICES Journal of Marine Science, 48: 281-290

  ! the data are given in Table 2 (p 286) and the four parameter model described on p287 is fitted:  note that the parameter 
  ! estimates differ from those on p288 because of rounding errors in the printing of the covariance matrices in Table 2


  IMPLICIT NONE

  INTEGER, PARAMETER    :: nhaul = 10

  INTEGER               :: i, j, k
  REAL (wp)             :: d(nnu,nnu), alpha(nalpha), se(nalpha)

  TYPE (HAUL_INFO)      :: haul(nhaul)


  ! read in parameter estimates and covariance matrices for each haul

  OPEN (FILE = 'nu.dat', UNIT = 9, STATUS = 'old')
  READ (9,*) ((haul(i)%nu(j), j = 1,nnu), ((haul(i)%r(j,k), k = 1,j), j = 1,nnu), i = 1, nhaul)


  ! construct x matrices

  DO i = 1, nhaul
    haul(i)%x = 0.0_wp
    haul(i)%x(1,1) = 1.0_wp
    haul(i)%x(2,2) = 1.0_wp
  END DO

  DO i = 7, nhaul
    haul(i)%x(1,3) = 1.0_wp
    haul(i)%x(2,4) = 1.0_wp
  END DO


  ! fit fixed and random effects selectivity model

  CALL FRSEL_FIT(haul, alpha = alpha, se = se, d = d)


  ! write output
 
  WRITE(*,100) (alpha(j), se(j), j = 1, nalpha), ((d(j,k), k = 1, 2), j = 1, 2)

  STOP

  100 FORMAT(/,'    alpha       se',/,4(/,2f9.3),//,'        d',/,2(/,2f9.4),//)

END PROGRAM EXAMPLE