A detailed derivation of improved ocean tidal equations in continuous (COTEs) and discrete (DOTEs) forms is presented. These equations feature the Boussinesq linear eddy dissipation law with a novel eddy viscosity that depends on the lateral mesh area, i.e., on mesh size and ocean depth. Analogously, the linear law of bottom friction is used with a new bottom friction coefficient depending on the bottom mesh area. The primary astronomical tide‐generating potential is modified by secondary effects due to the oceanic and terrestrial tides. The fully linearized equations are defined in a single‐layer ocean basin of realistic bathymetry varying from 50 m to 7,000 m. The DOTEs are set up on a 1oby 1ospherically graded grid system, using central finite differences in connection with Richardson's staggered computation scheme. Mixed single‐step finite differences in time are introduced, which enhance decay, dispersion, and stability properties of the DOTEs and facilitate—in Part II of this paper—a unique hydrodynamical interpolation of empirical tide data. The purely hydrodynamical modeling is completed by imposing boundary conditions consisting of no‐flow across and free‐slip along the mathematical ocean shorelines. Shortcomings of the constructed preliminary M2ocean tide charts are briefly discussed. Needed improvements of the model are left to Part II.