Merge reservoir 2 into reservoir 1
| Type | Intent | Optional | Attributes | Name | ||
|---|---|---|---|---|---|---|
| type(NeutralParticleReservoir), | intent(inout) | :: | res1 | Reservoir 1 | ||
| type(NeutralParticleReservoir), | intent(inout) | :: | res2 | Reservoir 2 | 
subroutine merge_reservoirs(res1,res2)
    !+ Merge reservoir 2 into reservoir 1
    type(NeutralParticleReservoir), intent(inout) :: res1
        !+ Reservoir 1
    type(NeutralParticleReservoir), intent(inout) :: res2
        !+ Reservoir 2
    integer :: i
    type(NeutralParticleReservoir) :: res3
    real(Float64) :: p, randomu(1)
    if(res1%n.le.0) then
        call init_reservoir(res1) !res2%k)
    endif
    if(res2%n.le.0) return
    !! If neither reservoir is full just fill res1 from res2 normally
    if ((res1%n.le.res1%k).and.(res2%n.le.res2%k)) then
        do i=1,res2%n
            call update_reservoir(res1, res2%R(i)%v, res2%R(i)%w)
        enddo
        res2 = res1
        return
    endif
    !! Switch res1 with res2 if res1 is not full but res2 is
    if (res1%n.le.res1%k) then
        res3 = res1
        res1 = res2
        res2 = res3
    endif
    !! if res2 is not full then fill res1 from res2 normally
    if (res2%n.le.res2%k) then
        do i=1,res2%n
            call update_reservoir(res1, res2%R(i)%v, res2%R(i)%w)
        enddo
    else !! if both are full select particles from res2 according to # particles seen
        p = float(res1%n)/(res1%n + res2%n)
        do i=1,res1%k
            call randu(randomu)
            if(randomu(1).gt.p) then
                res1%R(i) = res2%R(i)
            endif
        enddo
    endif
    !! set res1 to res2
    res2 = res1
end subroutine merge_reservoirs