Gradient-based algorithms used to solve aerodynamic shape optimization problems are regularly supported by the adjoint method (continuous or discrete) that can compute the gradients of the objective function with respect to the design variables. The continuous adjoint offers a physical insight into the adjoint equations and low memory footprint, but the selected discretization (usually based on intuition) affects the accuracy of the gradients. On the other hand, discrete adjoint computes sensitivity derivatives which are consistent with the discretized flow problem, with a higher memory footprint though. In this paper, the Think Discrete-Do Continuous (TDDC) adjoint is presented for turbulent and/or transitional flows. This bridges the gap between the two adjoint variants by proposing consistent discretization schemes for the continuous adjoint equations with a physical meaning; these are inspired by the hand-differentiated discrete adjoint. This work, extends a previous work of the group of authors, which was dealing with the Euler equations, by proposing consistent discretization schemes for turbulent and/or transitional flows. The development is performed on the in-house GPU-accelerated PUMA code. This uses a vertex-centered finite volume formulation on unstructured grids and second-order accurate discretization schemes (eg. Roe's scheme). The TDDC is verified in different grid sizes and applied for the shape optimization in turbulent/transitional external flows. For transitional cases, it is demonstrated that the impact of the ``frozen transition'' assumption can harm the accuracy of the computed gradients and should be avoided. The TDDC discretization schemes for the adjoint equations will be applied in industrial applications in a companion paper. This work is supported by the Hellenic Foundation for Research and Innovation (H.F.R.I.) under the “2nd Call for H.F.R.I. Research Projects to support Faculty Members & Researchers” (Project Number: 3821).