I proposed a change based on:
1) assume photoreduction is all UV
2) assume UV at surf is a constant fraction of PAR
(they diverge below the surface because UV is absorbed more readily by seawater)
3) get extinction coefficient from literature of UV that photoreduces DOM
4) from 2 and 3, compute penetration depth of UV dz_UV,
i.e. z where UV drops below some threshold
Because UV is absorbed so well, this will presumably be in the first model layer.
It will also be 0 when there is no light.
5) term that is added to DOM remin is dz_UV/dz(1)/(photoreduction timescale) * DOM

There are some tunable parameters here, but this approach has the following properties:
1) it turns off when there is no light
2) depth over which photoreduction occurs is independent of vertical grid.
3) as it is tied to UV absorption instead of PAR, it is focused in the very near surface, where UV is present.

Keith and Rob Letscher seemed to think this approach was better than what was in the code.
It still has ad-hoc assumptions. In particular, UV at high latitudes is probably not proportional to PAR.
I'm sure clouds attenuate UV and PAR differently too. But, it's something to work with.

Here's the math

UV(z=0) = C*PAR(z=0) ! C is unknown constant of proportionality
UV(z) = C*PAR(z=0)*exp(-z/L) ! L is unknown, but probably short, length scale
Let UV* denote threshold where photoreduction occurs only if UV>UV*.
(This is ad-hoc too. The process is probably a continuous function of UV, not a step function.)
Then z*, the depth at which photoreduction ceases (dz_UV from above), solves
UV* = C*PAR(z=0)*exp(-z*/L)
Solving for z* yields
z* = L*log(C*PAR(z=0)/UV*)
Note that only holds if z* is positive (i.e. argument of log()>1).
That is, if PAR(z=0) is not big enough, we don't want z* above the sea surface
So we really have
z* = L*log(C*PAR(z=0)/UV*), if PAR(z=0) > UV*/C
        0,     if PAR(z=0) < UV*/C
Now I really make stuff up. Other PAR based switches in the code use 1 W/m^2 as a threshold.
So I assume UV*/C = 1, out of convenience. So I get
z* = L*log(PAR(z=0)), if PAR(z=0) > 1
        0,     if PAR(z=0) < 1

Plugging this into my 5) from above yields that the additional remin from UV is
L*log(PAR(z=0))/dz(1)/(photoreduction timescale) * DOM (if PAR(z=0)>1)
I then redefine L to be L/10e2, which is dz(1) in our current grid and get
(L*log(PAR(z=0)))*(10.0e2/dz(1))/(photoreduction timescale) * DOM (if PAR(z=0)>1)
This way, 10e2/dz(1) is 1.0 for our current grid, but we get grid dependence built in.
Now I tune to get the same net photoreduction that Keith had previously tuned for.
I kept the photoreduction timescale to be the value that Keith had tuned for.
I did a diagnostic run where I computed the time space mean of log(max(PAR(z=0),1)).
I then solved for L so that L*mean(log(max(PAR(z=0),1)) = 1, to match Keith's values.
This yields L=0.4373.
Putting it all together yields the expression in the code:
(log(PAR_in(subcol_ind))*0.4373_r8) * (10.0e2_r8/dz1)