The complex-frequency-shifted perfectly matched layer (CFS-PML) technique can efficiently absorb near-grazing incident waves. In seismic wave modeling, CFS-PML has been implemented by the first-order-accuracy convolutional PML technique or second-order-accuracy recursive convolution PML technique. Both use different algorithms than the numerical scheme for the interior domain to update auxiliary memory variables in the PML and thus cannot be used directly with higher-order time-marching schemes. We work with an unsplit-field CFS-PML implementation using auxiliary differential equations (ADEs) to update the auxiliary memory variables. This ADE CFS-PML results in complete first-order differential equations. Thus, the numerical scheme for the interior domain can be used to solve ADE CFS-PML equations. We have implemented ADE CFS-PML in the finite-difference time-domain method and in anonstaggered-grid finite-difference method with the fourth-order Runge-Kutta scheme, demonstrating its straightforward implementation in different numerical time-marching schemes. We have also theoretically analyzed the role of the scalingfactor of CFS-PML; it transforms the PML to a transversely isotropic material, reducing the effective wave speed normal to the PML layer and bending the wavefront toward the normal direction of the PML layer. Our numerical tests indicate that the optimal value reduces the points per dominant wavelength at the outermost boundary to three, about half the value required by the numerical scheme. We also have found that the PML equations should be derived taking the free-surface boundary condition into account in finite-difference methods. Otherwise, the free surface in the PML layer causes instability or ineffective absorption of surface waves. Tests show that we can use a narrow-slice mesh with ADE CFS-PML to simulate full wave propagation efficiently in models with complex structure.