\textbf{The simple bird model}

New value code for alpha:

\begin{verbatim}
IF (alpha % buffer_index > SIZE(alpha % buffer) .OR. alpha % buffer_index < &
    1) THEN
    CALL generate('alpha', alpha % buffer)
    alpha % buffer_index = 1
END IF

alpha % one_dim(1) % new_value = alpha % buffer(alpha % buffer_index)
alpha % buffer_index = alpha % buffer_index + 1
\end{verbatim}

New values for functional parameters depending on alpha:

\begin{verbatim}
! calculating new value(s) for p

DO i2 = 1, 21
    DO j2 = 1, 4
        p % two_dim(i2, j2) % new_value = EXP(alpha % one_dim(1) % &
            new_value*northerness % one_dim(i2) % value)/(1+EXP(alpha % &
            one_dim(1) % new_value*northerness % one_dim(i2) % value))
    END DO
END DO
\end{verbatim}

Acceptation formula for alpha:

\begin{verbatim}
!  P(alpha') / P(alpha)

p_acc = 1
IF (alpha % one_dim(1) % new_value < -1.0_dp .OR. alpha % one_dim(1) % &
    new_value > 1.0_dp) THEN
    new_frequency = 0
ELSE
    new_frequency = 1.0 / (1.0_dp - (-1.0_dp))
END IF
IF (alpha % one_dim(1) % value < -1.0_dp .OR. alpha % one_dim(1) % value > &
    1.0_dp) THEN
    frequency = 0
ELSE
    frequency = 1.0 / (1.0_dp - (-1.0_dp))
END IF
p_acc = p_acc * new_frequency / frequency

! P( x | alpha') / P( x | alpha)

DO i=1, 21
    DO j=1, 4
        CALL user_bernoulli_freq(p % two_dim(i, j) % new_value, x % &
            two_dim(i, j) % value, new_frequency)
        CALL user_bernoulli_freq(p % two_dim(i, j) % value, x % two_dim(i, j) % &
            value, frequency)
        p_acc = p_acc * new_frequency / frequency
    END DO
END DO

! q( alpha', alpha) / q(alpha, alpha')

IF (alpha % one_dim(1) % value < -1.0_dp .OR. alpha % one_dim(1) % value > &
    1.0_dp) THEN
    frequency = 0
ELSE
    frequency = 1.0 / (1.0_dp - (-1.0_dp))
END IF
IF (alpha % one_dim(1) % new_value < -1.0_dp .OR. alpha % one_dim(1) % &
    new_value > 1.0_dp) THEN
    new_frequency = 0
ELSE
    new_frequency = 1.0 / (1.0_dp - (-1.0_dp))
END IF
p_acc = p_acc * frequency / new_frequency
\end{verbatim}

Acceptation code for alpha (if the update strategy is sequential): 

\begin{verbatim}
IF (p_acc > prob(prob_index)) THEN
    alpha % one_dim(1) % value = alpha % one_dim(1) % new_value
    alpha % one_dim(1) % successful_changes = alpha % one_dim(1) % &
        successful_changes + 1
    
    DO i2 = 1, 21
        DO j2 = 1, 4
            p % two_dim(i2, j2) % value = p % two_dim(i2, j2) % new_value
        END DO
    END DO
END IF
prob_index = prob_index + 1
IF (prob_index > SIZE(prob) .OR. prob_index < 1) THEN
    CALL G05FAF(0.0_dp, 1.0_dp, SIZE(prob), prob)
    prob_index = 1
END IF
\end{verbatim}

Acceptation code for alpha (if the update strategy is random): 

\begin{verbatim}
IF (p_acc > prob(prob_index)) THEN
    alpha % one_dim(1) % value = alpha % one_dim(1) % new_value
    alpha % one_dim(1) % successful_changes = alpha % one_dim(1) % &
        successful_changes + 1
    
    IF (alpha % one_dim(1) % successful_changes == alpha % one_dim(1) % &
        update_count) THEN
        parameters_achieved = parameters_achieved + 1
    END IF
    DO i2 = 1, 21
        DO j2 = 1, 4
            p % two_dim(i2, j2) % value = p % two_dim(i2, j2) % new_value
        END DO
    END DO
END IF
prob_index = prob_index + 1
IF (prob_index > SIZE(prob) .OR. prob_index < 1) THEN
    CALL G05FAF(0.0_dp, 1.0_dp, SIZE(prob), prob)
    prob_index = 1
END IF
\end{verbatim}

New value code for x:

\begin{verbatim}
IF (x % buffer_index > SIZE(x % buffer) .OR. x % buffer_index < 1) THEN
    CALL generate('x', x % buffer)
    x % buffer_index = 1
END IF

x % two_dim(i, j) % new_value = x % buffer(x % buffer_index)
x % buffer_index = x % buffer_index + 1
\end{verbatim}

New values for functional parameters depending on x:

Nothing should be printed.

Acceptation formula for x:

\begin{verbatim}
!  P(x') / P(x)

p_acc = 1
CALL user_bernoulli_freq(p % two_dim(i, j) % value, x % two_dim(i, j) % &
    new_value, new_frequency)
CALL user_bernoulli_freq(p % two_dim(i, j) % value, x % two_dim(i, j) % &
    value, frequency)
p_acc = p_acc * new_frequency / frequency

! P( obs | x') / P( obs | x)

CALL user_defined_points_freq(x % two_dim(i, j) % new_value, obs % two_dim(i, &
    j) % value, new_frequency)
CALL user_defined_points_freq(x % two_dim(i, j) % value, obs % two_dim(i, j) &
    % value, frequency)
p_acc = p_acc * new_frequency / frequency

! q( x', x) / q(x, x')

IF (x % two_dim(i, j) % value < 0 .OR. x % two_dim(i, j) % value > 1) THEN
    frequency = 0
ELSE
    frequency = 1.0 / (REAL(1) - (REAL(0)) + 1.0)
END IF
IF (x % two_dim(i, j) % new_value < 0 .OR. x % two_dim(i, j) % new_value > 1) &
    THEN
    new_frequency = 0
ELSE
    new_frequency = 1.0 / (REAL(1) - (REAL(0)) + 1.0)
END IF
p_acc = p_acc * frequency / new_frequency
\end{verbatim}

Acceptation code for x (if the update strategy is sequential): 

\begin{verbatim}
IF (p_acc > prob(prob_index)) THEN
    x % two_dim(i, j) % value = x % two_dim(i, j) % new_value
    x % two_dim(i, j) % successful_changes = x % two_dim(i, j) % &
        successful_changes + 1
    
END IF
prob_index = prob_index + 1
IF (prob_index > SIZE(prob) .OR. prob_index < 1) THEN
    CALL G05FAF(0.0_dp, 1.0_dp, SIZE(prob), prob)
    prob_index = 1
END IF
\end{verbatim}

Acceptation code for x (if the update strategy is random): 

\begin{verbatim}
IF (p_acc > prob(prob_index)) THEN
    x % two_dim(i, j) % value = x % two_dim(i, j) % new_value
    x % two_dim(i, j) % successful_changes = x % two_dim(i, j) % &
        successful_changes + 1
    
    IF (x % two_dim(i, j) % successful_changes == x % two_dim(i, j) % &
        update_count) THEN
        parameters_achieved = parameters_achieved + 1
    END IF
END IF
prob_index = prob_index + 1
IF (prob_index > SIZE(prob) .OR. prob_index < 1) THEN
    CALL G05FAF(0.0_dp, 1.0_dp, SIZE(prob), prob)
    prob_index = 1
END IF
\end{verbatim}

\textbf{The simple bird model with spatial dependence}

Only sections that differ from the previous test case are described here.

New values for functional parameters depending on alpha:

\begin{verbatim}
! calculating new value(s) for p

DO i2 = 1, 21
    DO j2 = 1, 4
        p % two_dim(i2, j2) % new_value = EXP(alpha % one_dim(1) % &
            new_value*northerness % one_dim(i2) % value+q % two_dim(i2, j2) % &
            value)/(1+EXP(alpha % one_dim(1) % new_value*northerness % & 
            one_dim(i2) % value+q % two_dim(i2, j2) % value))
    END DO
END DO
\end{verbatim}

Acceptation code for x (if the update strategy is sequential): 

\begin{verbatim}
IF (p_acc > prob(prob_index)) THEN
    x % two_dim(i, j) % value = x % two_dim(i, j) % new_value
    x % two_dim(i, j) % successful_changes = x % two_dim(i, j) % &
        successful_changes + 1
    
    ! 1 to number of neighbours
    DO k = 1, spatial(i, 1)
        sum_int = 0
        neighbour = spatial(i, k + 1)
        ! updating q % two_dim(spatial(i, k), j)
        DO h = 1, spatial(neighbour, 1)
            sum_int = sum_int + x % two_dim(spatial(neighbour, h + 1), j) % &
                value
        END DO
        q % two_dim(neighbour, j) % value = sum_int
    END DO
    DO k = 1, spatial(i, 1)
        neighbour = spatial(i,  k +1)
        p % two_dim(neighbour, j) % value = EXP(alpha % one_dim(1) % &
            value*northerness % one_dim(neighbour) % value+q % &
            two_dim(neighbour, j) % value)/(1+EXP(alpha % one_dim(1) % &
            value*northerness % one_dim(neighbour) % value+q % &
            two_dim(neighbour, j) % value))
    END DO
END IF
prob_index = prob_index + 1
IF (prob_index > SIZE(prob) .OR. prob_index < 1) THEN
    CALL G05FAF(0.0_dp, 1.0_dp, SIZE(prob), prob)
    prob_index = 1
END IF
\end{verbatim}

Acceptation code for x (if the update strategy is random): 

\begin{verbatim}
IF (p_acc > prob(prob_index)) THEN
    x % two_dim(i, j) % value = x % two_dim(i, j) % new_value
    
    x % two_dim(i, j) % successful_changes = x % two_dim(i, j) % &
        successful_changes + 1
    
    IF (x % two_dim(i, j) % successful_changes == x % two_dim(i, j) % &
        update_count) THEN
        parameters_achieved = parameters_achieved + 1
    END IF
      
    ! 1 to number of neighbours
    DO k = 1, spatial(i, 1)
        sum_int = 0
        neighbour = spatial(i, k + 1)
        ! updating q % two_dim(spatial(i, k), j)
        DO h = 1, spatial(neighbour, 1)
            sum_int = sum_int + x % two_dim(spatial(neighbour, h + 1), j) % &
                value
        END DO
        q % two_dim(neighbour, j) % value = sum_int
    END DO
    DO k = 1, spatial(i, 1)
        neighbour = spatial(i,  k +1)
        p % two_dim(neighbour, j) % value = EXP(alpha % one_dim(1) % &
            value*northerness % one_dim(neighbour) % value+q % &
            two_dim(neighbour, j) % value)/(1+EXP(alpha % one_dim(1) % &
            value*northerness % one_dim(neighbour) % value+q % &
            two_dim(neighbour, j) % value))
    END DO
END IF
prob_index = prob_index + 1
IF (prob_index > SIZE(prob) .OR. prob_index < 1) THEN
    CALL G05FAF(0.0_dp, 1.0_dp, SIZE(prob), prob)
    prob_index = 1
END IF
\end{verbatim}