Code:

!
! ARHN_MLE_NEG_LOGLF
!
subroutine arhn_mle_neg_loglf(num_parm,mle_theta,neg_loglf)
! Define the output variables
double precision, intent(out) :: neg_loglf
! Define the input variables
integer, intent(in) :: num_parm
double precision, intent(in) :: mle_theta(num_parm)
! Define the local variables
double precision :: loglf
double precision :: loc_alpha(0:2)
double precision :: loc_beta
double precision :: loc_theta(1:num_ntot)
! Calculate the local variables
loc_alpha(0:2) = mle_theta(1:3)
loc_beta = mle_theta(4)
loc_theta = exp(loc_alpha(0) - loc_alpha(1) * vect_t * exp(loc_alpha(2) * vect_v))
! Calculate the log likelihood function
loglf = (sum(n_dp) * log(loc_beta)) - (loc_beta * sum(loc_alpha(0) - loc_alpha(1) * vect_t * exp(loc_alpha(2) * vect_v))) + &
& ((loc_beta - 1d0) * sum(log(vect_x))) - (sum((vect_x / loc_theta)**loc_beta))
! Calculate the negative log likelihood function
neg_loglf = -loglf
endsubroutine arhn_mle_neg_loglf
!
! ARHN_MLE_NEG_LOGLF_1DER
!
subroutine arhn_mle_neg_loglf_1der(num_parm,mle_theta,neg_loglf_1der)
! Define the output variables
double precision, intent(out) :: neg_loglf_1der(num_parm)
! Define the input variables
integer, intent(in) :: num_parm
double precision, intent(in) :: mle_theta(num_parm)
! Define the local variables
double precision :: loglf_1der(num_parm)
double precision :: loc_alpha(0:2)
double precision :: loc_beta
double precision :: loc_theta(1:num_ntot)
! Calculate the local variables
loc_alpha(0:2) = mle_theta(1:3)
loc_beta = mle_theta(4)
loc_theta = exp(loc_alpha(0) - loc_alpha(1) * vect_t * exp(loc_alpha(2) * vect_v))
! Calculate the first order partial derivatives of the log likelihood function
loglf_1der(1) = -loc_beta * sum(n_dp) + loc_beta * sum((vect_x / loc_theta)**loc_beta)
loglf_1der(2) = loc_beta * sum(vect_t * exp(loc_alpha(2) * vect_v) * (1d0 - (vect_x / loc_theta)**loc_beta))
loglf_1der(3) = loc_beta * sum(vect_t * exp(loc_alpha(2) * vect_v) * (1d0 - (vect_x / loc_theta)**loc_beta) * loc_alpha(1) * vect_v)
loglf_1der(4) = sum((1d0 / loc_beta) + log(vect_x / loc_theta) * (1d0 - (vect_x / loc_theta)**loc_beta))
! Calculate the first order partial derivatives of the negative log likelihood function
neg_loglf_1der = -loglf_1der
endsubroutine arhn_mle_neg_loglf_1der
!
! ARHN_MLE_NEG_LOGLF_2DER
!
subroutine arhn_mle_neg_loglf_2der(num_parm,mle_theta,neg_loglf_2der,LDH)
! Define the output variables
double precision, intent(out) :: neg_loglf_2der(LDH,num_parm)
! Define the input variables
integer, intent(in) :: LDH
integer, intent(in) :: num_parm
double precision, intent(in) :: mle_theta(num_parm)
! Define the local variables
double precision :: loglf_2der(LDH,num_parm)
double precision :: loc_alpha(0:2)
double precision :: loc_beta
double precision :: loc_theta(1:num_ntot)
! Calculate the local variables
loc_alpha(0:2) = mle_theta(1:3)
loc_beta = mle_theta(4)
loc_theta = exp(loc_alpha(0) - loc_alpha(1) * vect_t * exp(loc_alpha(2) * vect_v))
! Calculate the second order partial derivatives of the log likelihood function
loglf_2der(1,1) = -(loc_beta**2) * sum((vect_x / loc_theta)**loc_beta)
loglf_2der(1,2) = (loc_beta**2) * sum(vect_t * exp(loc_alpha(2) * vect_v) * (vect_x / loc_theta)**loc_beta)
loglf_2der(1,3) = (loc_beta**2) * sum(vect_t * exp(loc_alpha(2) * vect_v) * loc_alpha(1) * vect_v * (vect_x / loc_theta)**loc_beta)
loglf_2der(1,4) = -sum(n_dp) + sum((1d0 + loc_beta * log(vect_x / loc_theta)) * (vect_x / loc_theta)**loc_beta)
loglf_2der(2,1) = loglf_2der(1,2)
loglf_2der(2,2) = -(loc_beta**2) * sum((vect_t * exp(loc_alpha(2) * vect_v))**2 * (vect_x / loc_theta)**loc_beta)
loglf_2der(2,3) = -(loc_beta**2) * sum((vect_t * exp(loc_alpha(2) * vect_v))**2 * loc_alpha(1) * vect_v * (vect_x / loc_theta)**loc_beta) + loc_beta * sum(vect_t * vect_v * exp(loc_alpha(2) * vect_v) * (1d0 - (vect_x / loc_theta)**loc_beta))
loglf_2der(2,4) = sum((vect_t * exp(loc_alpha(2) * vect_v)) * (1d0 - (1d0 + loc_beta * log(vect_x / loc_theta)) * (vect_x / loc_theta)**loc_beta))
loglf_2der(3,1) = loglf_2der(1,3)
loglf_2der(3,2) = loglf_2der(2,3)
loglf_2der(3,3) = -(loc_beta**2) * sum((vect_t * exp(loc_alpha(2) * vect_v) * loc_alpha(1) * vect_v)**2 * (vect_x / loc_theta)**loc_beta) + loc_beta * sum(vect_t * exp(loc_alpha(2) * vect_v) * (1d0 - (vect_x / loc_theta)**loc_beta) * loc_alpha(1) * vect_v**2)
loglf_2der(3,4) = sum((vect_t * exp(loc_alpha(2) * vect_v) * loc_alpha(1) * vect_v) * (1d0 - (1d0 + loc_beta * log(vect_x / loc_theta)) * (vect_x / loc_theta)**loc_beta))
loglf_2der(4,1) = loglf_2der(1,4)
loglf_2der(4,2) = loglf_2der(2,4)
loglf_2der(4,3) = loglf_2der(3,4)
loglf_2der(4,4) = -sum(n_dp) / (loc_beta**2) - sum((log(vect_x / loc_theta))**2 * (vect_x / loc_theta)**loc_beta)
! Calculate the second order partial derivatives of the negative log likelihood function
neg_loglf_2der = -loglf_2der
endsubroutine arhn_mle_neg_loglf_2der

The global variables come up because n_dp, vect_x, vect_t, and vect_v are defined in the main program, but not the internal subroutines. I understand that the manual says the subroutines need to be EXTERNAL. I tried several different ways, but I couldn't determine any efficient method of passing the data to the subroutine unless I used an internal subroutine with global variables. In my test programs (all based on the sample problem in the reference manual), the results were the same for my method.