& rhoConstFresh, rho1Ref, tRef, sRef, phiRef, dBdrRef,

      _RL phiRef(2*Nr+1)

        hRef(K) = phiRef(K)*recip_gravity

          deepFacC(k)  = ( rSphere + phiRef(2*k)*recip_gravity

          deepFacF(k)  = ( rSphere + phiRef(2*k-1)*recip_gravity

       phiRef(1) = 0. _d 0

          phiRef(k+1) = phiRef(k)+ddPI*tRef(ks)

          phiRef(2*k)   = phiRef(1) + ddPI*tRef(k)

          phiRef(2*k+1) = phiRef(2*k) + ddPI*0.5*tRef(k)

          phiRef(2*k+2) = phiRef(2*k)

          phiRef(2*k+1) = phiRef(2*k) + ddPI*tRef(k)

        phiRef(1) = 0. _d 0

          phiRef(2*k)   = gravitySign*gravity*(rC(k) - Ro_SeaLevel)

          phiRef(2*k+1) = gravitySign*gravity*(rF(k+1)-Ro_SeaLevel)

        phiRef(1) = -recip_rhoConst*rF(1)

          phiRef(2*k)   = -recip_rhoConst*rC(k)

          phiRef(2*k+1) = -recip_rhoConst*rF(k+1)

     &    phiRef(k)*recip_gravity

        phiRef(k) = 0.

          tmpVar(k) = phiRef(2*k-1)

          tmpVar(k) = phiRef(2*k)