Approximate diffs after 5-day runs (bug is in oplus.F):
2011 day 160: Weimer potential model, IMF_OMNI data, GPI run

The original bug in oplus.F was that it used the wrong reduced
temperature (tpj=(te+ti)) instead of tr=(tn+ti)) to calculate 
the ambi-polar diffusion (sub rrk). After this the bugfix and
nobugfix codes are effectively the same.

Field TPTR_rrk in bugfix:   tr  = 0.5*(tn+ti) (correct reduced temperature )
Field TPTR_rrk in nobugfix: tpj = 0.5*(te+ti) (wrong)

%Diffs = (bugfix-nobugfix)/nobugfix*100

Field		%Diffs
-----------------------------------------------
TPTR_rrk	-1372%    tr(tn+ti) - tpj(te+ti)
                          bugfix      nobugfix

DJ		40%          ambipolar diffusion
DIFFJ		-14->8%      output of diffus
TPJ		-10->20%     plasma temp (te+ti) * O+
TP1		-3.2->4.6%   0.5*(te+ti)
XIOP2D		-10->15%     sources or sinks (also XIOP2P)
DJINT		40%          DJ on interfaces
OP_LOSS		-7->7%       O+ loss rate
OPOUT		-12->18%     O+ output
NE		-14->12% 

----------------------------
Old Code (nobugfix/oplus.F):

tpj = 0.5*(te+ti)        ! incorrect (should be tr=0.5*(tn+ti))
call rrk(...,tpj,...)    ! ambipolar diffusion output DJ

tpj = 2.*tpj             ! this is now the correct plasma temp (te+ti) 
call diffus(...,tpj,...) ! DIFFJ output

tpj = tpj * op           ! TPJ diag
call bdotdh(...,tpj,...) ! BDOTDH output

tp = 0.5*(ti+te)         ! TP1 diag: now using average rather than total?
tphdz0,1 = ...tp...      ! to be used later for p,q,r coeffs

DJINT = (DJ on interfaces) ! to be used later for p,q,r coeffs

XIOP2D,XIOP2P sources and sinks
OP_LOSS

trsolv()
filter_op()

OPOUT

----------------------------
New Code (bugfix/oplus.F):

tr = 0.5*(tn+ti)         ! reduced temperature
call rrk(...,tr,...)     ! ambipolar diffusion output DJ

tp = (te+ti)             ! plasma temperature
call diffus(...,tp,...)  ! DIFFJ output

tp = tp * op             ! TPJ diag
call bdotdh(...,tp,...)  ! BDOTDH output

tp1 = 0.5*(ti+te)        ! TP1 diag: now using average rather than total?
tphdz0,1 = ...tp1...     ! to be used later for p,q,r coeffs

DJINT = (DJ on interfaces) ! to be used later for p,q,r coeffs

XIOP2D,XIOP2P sources and sinks
OP_LOSS

trsolv()
filter_op()
enforce_opfloor...

OPOUT

----------------------------