O rhoP0,
& - rhoP0(i,j)*p1*dKdtheta )
_RL rhoP0 (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
& - rhoP0(i,j)*p1*dKdS )
rhoP0 = rfresh + rsalt
rhoLoc = rhoP0/(1. _d 0 - p1/bulkMod)
rhoLoc(i,j) = rhoP0(i,j)
_RL rhoP0(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
rhoP0(i,j) = rfresh + rsalt
rhoP0(i,j) = 0. _d 0
_RL rfresh, rsalt, rhoP0
rhoP0 = 0. _d 0