Return union and complement of three sets I would be especially happy to hear opinions about complement_union. 
The rest of the code is just there to test the output.

program test_complement_union
    implicit none(type, external)

    call assert(all([3, 5] == complement_union([1, 2, 3], [1, 2], [5])))
    call assert(all([2, 3, 4] == complement_union([1, 2, 3], [1], [4])))
    call assert(all([1, 3, 5] == complement_union([1, 2, 3], [2], [5])))
contains
    !&gt; Merge C into A and remove values that are in B.
    !&gt; The result can be written with set notation as A ∪ C / B.
    !&gt; Preconditions (not tested!):
    !&gt;      1. B is a subset of A
    !&gt;      2. A and C are disjoint
    !&gt;      3. B and C are disjoint
    !&gt;      4. A, B, and C are sorted.
    pure function complement_union(A, B, C) result(D)
        integer, intent(in) :: A(:), B(:), C(:)
        integer :: D(size(A) + size(C) - size(B))

        ! We will index as A(i), B(j), C(k), D(l)
        integer :: i, j, k, l

        i = 1
        j = 1
        k = 1
        l = 1
        do while(l &lt;= size(D))
            ! Only indices from C have to be added to A
            ! We use assumption that B is a subset of A
            if (i &gt; size(A)) then
                D(l) = C(k)
                k = k + 1
                l = l + 1
            ! Neither indices from B have to be deleted in A
            ! nor indices from C have to be added from C to A.
            else if (j &gt; size(B) .and. k &gt; size(C)) then
                D(l) = A(i)
                i = i + 1
                l = l + 1
            ! No more indices from B have to be deleted in A
            else if (j &gt; size(B)) then
                if (A(i) &lt; C(k)) then
                    D(l) = A(i)
                    i = i + 1
                    l = l + 1
                else if (A(i) &gt; C(k)) then
                    D(l) = C(k)
                    k = k + 1
                    l = l + 1
                end if
            ! No more indices have to be added from C to A
            else if (k &gt; size(C)) then
                if (A(i) /= B(j)) then
                    D(l) = A(i)
                    i = i + 1
                    l = l + 1
                else
                    i = i + 1
                    j = j + 1
                end if
            ! Normal case:
            ! Merge C sorted into A excluding values from B.
            else if (A(i) &lt; C(k)) then
                if (A(i) /= B(j)) then
                    D(l) = A(i)
                    i = i + 1
                    l = l + 1
                else
                    i = i + 1
                    j = j + 1
                end if
            else if (A(i) &gt; C(k)) then
                if (A(i) /= B(j)) then
                    D(l) = C(k)
                    k = k + 1
                    l = l + 1
                else
                    D(l) = C(k)
                    i = i + 1
                    j = j + 1
                    k = k + 1
                    l = l + 1
                end if
            end if
        end do
    end function

    pure subroutine assert(cond)
        logical, intent(in) :: cond
        if (.not. cond) error stop
    end subroutine
end program
```
This code has two limitations, as I see it:

in treating in a single go the combination of three sets, it gives rise to a high cyclomatic complexity: a dozen if statements (including the while) in a single function. This circumstance limits the readability and maintainability of the code. It also makes it less extensible to more general cases.
the assumption that B is a subset of A and C is disjoint from either is very particular, and likely to limit the reuse of the code elsewhere.

As a minor point, the assertions stop the program in case a problem is detected but they do not indicate which test fails, or if other tests would have failed, which comes in handy in debug.
The input vectors are assumed to be sorted in increasing order. If that were not the case, one could do that first, and increase the generality of the code further.
Leaving aside for a moment the possible generalization of the code, by splitting the operation into a subtraction of B from A and a subsequent addition of C, the complexity drops drastically. Instead of a single pure function with 55 lines of code (header and declarations excluded) and 12 branching points (if and while), one can write three pure functions, with 2, 11, and 21 lines of code and only 0, 2, and 3 branching points, respectively:
program test_complement_union
  
    implicit none(type, external)
    
    call assert( &quot;t1&quot;, all([3, 5]       == complement_union([1, 2, 3], [1, 2], [5]      )))
    call assert( &quot;t2&quot;, all([2, 3, 4]    == complement_union([1, 2, 3], [1]   , [4]      )))
    call assert( &quot;t3&quot;, all([1, 3, 5]    == complement_union([1, 2, 3], [2]   , [5]      )))
    call assert( &quot;t4&quot;, all([1, 4, 5, 6] == complement_union([1, 2, 3], [2, 3], [4, 5, 6])))

  contains

    !&gt; Merge C into A and remove values that are in B.
    !&gt; The result can be written with set notation as A ∪ C / B.
    !&gt; Preconditions (not tested!):
    !&gt;      1. B is a subset of A
    !&gt;      2. A and C are disjoint
    !&gt;      3. B and C are disjoint
    !&gt;      4. A, B, and C are sorted.
    pure function complement_union(A, B, C) result(D)
      integer, intent(in) :: A(:), B(:), C(:)
      integer :: A_B( size(A) - size(B) )
      integer :: D  ( size(A) + size(C) -size(B) )
      A_B  = subtract_subset(A,B)
      D    = union_disjoint_sets(A_B,C)
    end function complement_union
    
    !&gt; Subtract a proper ordered subset from a given ordered set
    pure function subtract_subset( A, B ) result( C )
      integer, intent(in) :: A(:), B(:)
      integer :: C ( size(A) - size(B) ), iA, iB, iC
      iC = 0
      iA = 1
      !.. All elements of B are expected to be found in A
      do iB = 1, size(B)
         !.. Copies the elements of A between consecutive
         !   occurrences of B elements
         do while( A(iA) /= B(iB) )
            iC = iC + 1
            C(iC) = A(iA)
            iA = iA + 1
         enddo
         !.. Skips the B element
         iA = iA + 1
      enddo 
      !.. Copies the elements of A past the last element of B
      if( iA &lt;= size(A) ) C(iC+1:) = A(iA:)
    end function subtract_subset
    
    !&gt; Computes the union of non-empty disjoint ordered sets
    pure function union_disjoint_sets( A, B ) result( C )
      integer, intent(in) :: A(:), B(:)
      integer :: C( size(A) + size(B) )
      integer :: iA, iB, iC, nA, nB
      nA = size(A)
      nB = size(B)
      iA = 1
      iB = 1
      iC = 0
      C  = 0
      !.. Takes an element from either A or B
      !   until one of the sets is empty
      do while( iA &lt;= nA .and. iB &lt;= nB )
         iC = iC + 1
         if( A(iA) &lt;= B(iB) )then
            C(iC) = A(iA)
            iA = iA + 1
         else
            C(iC) = B(iB)
            iB = iB + 1
         endif
      end do
      !.. Copies to C the residual elements in
      !   the remaining non-empty set
      if( iA &lt;= nA )then
         C(iC+1:)=A(iA:)
      else
         C(iC+1:)=B(iB:)
      endif
    end function union_disjoint_sets
    
    subroutine assert(msg,cond)
      use, intrinsic :: iso_fortran_env, only : OUTPUT_UNIT
      character(len=*), intent(in) :: msg
      logical         , intent(in) :: cond
      write(OUTPUT_UNIT,&quot;(a)&quot;,advance=&quot;no&quot;) msg//&quot; - &quot;
      if( cond )then
         write(OUTPUT_UNIT,&quot;(a)&quot;) &quot;passed&quot;
      else
         write(OUTPUT_UNIT,&quot;(a)&quot;) &quot;FAILED&quot;
      endif
    end subroutine assert

  end program test_complement_union

If the restrictions on the three sets, A, B, and C are removed, except for their being ordered, in general, the size of ( A U C ) / B will be different from size(A)-size(B)+size(C). Furthermore, ( A / B ) U C /= ( A U C ) / B. Let's consider the latter operation, which is a bit more challenging. If one does not want give up automatic arrays, then I think there is no alternative to having trailing &quot;empty&quot; entries in the result (e.g., zeros). Even in this case, it is best to split the work into functions that provide two disjoint sets to be merged. The subtraction here is more complicated because in general B would not be a subset of A and C may have elements in common with both B and A. Therefore, one must first figure out the A/B, C/B and C/B/A, by means of a difference function, and then merge the result, which is simple, since the sets are disjoint :
$$( A \cup C ) / B = ( A / B ) \cup ( ( C / B ) / A ) $$
where $$( A / B ) \cap ( ( C / B ) / A ) = \varnothing$$
Here's a sketch of the idea
program test_complement_union
  
    implicit none(type, external)
    
    call assert( &quot;t1&quot;, all([3, 5, 0, 0]       == complement_union([1, 2, 3], [1, 2], [5]      )))
    call assert( &quot;t2&quot;, all([2, 3, 4, 0]       == complement_union([1, 2, 3], [1]   , [4]      )))
    call assert( &quot;t3&quot;, all([1, 3, 5, 0]       == complement_union([1, 2, 3], [2]   , [5]      )))
    call assert( &quot;t4&quot;, all([1, 3, 5, 0, 0, 0] == complement_union([1, 2, 3], [2, 9], [1, 2, 5])))

  contains

    !.. Assumes input vectors are sorted (!)
    !   If that is not the case, the vectors should
    !   be sorted into auxiliary arrays
    !..
    function complement_union(A, B, C) result(D)
      integer, intent(in) :: A(:), B(:), C(:)
      integer :: A_B(size(A))
      integer :: C_B(size(C))
      integer :: C_AB(size(C))
      integer :: D (size(A) + size(C))
      A_B  = difference(A,B)
      C_B  = difference(C,B)
      C_AB = difference(C_B,A_B)
      D    = union_disjoint_sets(A_B,C_AB)
    end function complement_union
    
    !.. Computes the difference between the two sets,
    !   generating a set with the same size as the first
    !   with trailing zeros for the erased fields.
    function difference( A, B ) result( C )
      integer, intent(in) :: A(:), B(:)
      integer :: C ( size(A) ), iA, iB, iC, nA, nB, nC
      nA=count(A/=0)
      nB=count(B/=0)
      if(nB==0)then
         C=A
         return
      endif
      iC=0
      iB=1
      nC=nA
      outer: do iA = 1, nA
         do while( B(iB) &lt; A(iA) )
            iB = iB + 1
            if( iB &gt; nB) exit outer
         end do
         if( B(iB) == A(iA) )then
            nC=nC-1
            cycle outer
         endif
         iC = iC + 1
         C(iC) = A(iA)
      enddo outer
      if( iB &gt; nB ) C(iC+1:nC) = A(iA:nA)
      if( nC &lt; size(A) ) C(nC+1:) = 0
    end function difference
    
    !.. Computes the union of disjoint sets
    pure function union_disjoint_sets( A, B ) result( C )
      integer, intent(in) :: A(:), B(:)
      integer :: C( size(A) + size(B) )
      integer :: iA, iB, iC, nA, nB, nC
      nA = count( A /= 0 )
      nB = count( B /= 0 )
      nC = nA + nB
      iA = 1
      iB = 1
      iC = 0
      C  = 0
      do while( iA &lt;= nA .and. iB &lt;= nB )
         iC = iC + 1
         if( A(iA) &lt;= B(iB) )then
            C(iC) = A(iA)
            iA = iA + 1
         else
            C(iC) = B(iB)
            iB = iB + 1
         endif
      end do
      if( iA &lt;= nA )then
         C(iC+1:nC)=A(iA:nA)
      else
         C(iC+1:nC)=B(iB:nB)
      endif
    end function union_disjoint_sets
    
    subroutine assert(msg,cond)
      use, intrinsic :: iso_fortran_env, only : OUTPUT_UNIT
      character(len=*), intent(in) :: msg
      logical         , intent(in) :: cond
      write(OUTPUT_UNIT,&quot;(a)&quot;,advance=&quot;no&quot;) msg//&quot; - &quot;
      if( cond )then
         write(OUTPUT_UNIT,&quot;(a)&quot;) &quot;passed&quot;
      else
         write(OUTPUT_UNIT,&quot;(a)&quot;) &quot;FAILED&quot;
      endif
    end subroutine assert

  end program test_complement_union

PS:
For folks not familiar with the simple syntax of modern Fortran, I suggest the following quick modern Fortran tutorial.