if I want to display the current value (pre-update), I run this code:
        print *, U(:,i,j), 'before'
and get the expected output:
   4.000000       4.000000      0.0000000E+00   3.785714     before
the result iterates through a range of i's and j's and spits out many
reasonable numbers

When I want to display the current value along with what the updated
solution should be, I run this code:
        print *, U(:,i,j), 'before'
        print *,  U(:,i,j)-dt*jac(i,j)*((Fp-Fm)+(Gp-Gm)), 'update'
and get this (expected) output, again, only one item in a long list of
numbers is pasted:
   4.000000       4.000000      0.0000000E+00   3.785714     before
   4.004312       4.004308      0.0000000E+00   6.184789     update

Now, if I try to reassign the vector with the new solution, I use this
code:
        print *, U(:,i,j), 'before'
        print *,  U(:,i,j)-dt*jac(i,j)*((Fp-Fm)+(Gp-Gm)), 'update'
        U(:,i,j)= U(:,i,j)-dt*jac(i,j)*((Fp-Fm)+(Gp-Gm))

I get this output:
   4.000000       4.000000      0.0000000E+00   3.785714     before
 NaN            NaN            NaN            NaN            update

The updated value got changed to NaN.  This does not happen for all
iterations, some are NaN, others are normal-looking values, and it
appears seemingly at random.  I even made an explicit loop to do the
reassigning, like so:

        do k=1,4
          U(k,i,j)= U(k,i,j) - dt*jac(i,j)*((Fp(k)-Fm(k))
     &              +(Gp(k)-Gm(k)))

and the output is the same.  Do you have any idea what I am doing
incorrectly?  I apologize for not being able to paste my whole code...
its length does not permit it.

Exceeding array bounds can also do all kinds of strange things. As you
do have arrays, and no evidence relating to teh validity of the index
values, that possibility can't be ignored.