saltFreeze = (- bqe - SQRT(discrim))*recip_aqe
IF ( saltFreeze .LT. 0. _d 0 )
& saltFreeze = (- bqe + SQRT(discrim))*recip_aqe
thetaFreeze = a0*saltFreeze + eps4
_RL thetaFreeze, saltFreeze
& * ( 1. _d 0 - sLoc(I,J) / saltFreeze )
& * ( saltFreeze - sLoc(I,J) )
_RL thetaFreeze, saltFreeze, recip_Cp