This is the function I use, I've added some comments:
FUNCTION ModelineCreate ( BYVAL xres AS LONG, BYVAL yres AS LONG, BYVAL vfreq AS SINGLE, BYVAL rotation AS LONG, index AS LONG ) AS LONG
  LOCAL i, j, result, interlace, margen, DotClockIdx AS LONG
  LOCAL hhh, hhi, hhf, hht, newHhh, newHhi, newHhf, newHht, vvv, vvi, vvf, vvt, vIncr, vfreqLabel AS LONG
  LOCAL hfreq, Dotclock, DotclockReq, vvtIni, VBlankLines, diff, newDiff, newVfreq, hfreqReal, vfreqReal AS SINGLE
' We set default 'interlace' to 1 (no interlace). When set to 2 it means interlace is on. I use these values as they are
' useful for calculations, instead of 0, 1.
  INCR interlace
' This part is just to filter max and min resolutions, hfreq and vfreq values to keep them inside our monitor limits.
' It also handles vertical games (rotated resolutions). But its focused on lowres monitors with some hardcoded conditionals,
' so you may want to modify things to adapt it to multifrequency monitors. However, this should better be kept outside
' this function.
  IF xres < XresMin THEN xres = XresMin : result = result OR %SimilarRes
  IF yres < YresMin THEN yres = YResMin : result = result OR %SimilarRes
  IF rotation = 1 THEN VerticalToHorizontal ( xres, yres )
  IF vfreq < VfreqMin THEN
    IF vfreq * 2 <= VfreqMax THEN
      vfreq = vfreq * 2
      result = result OR %DoubleVfreq
    ELSE
      vfreq = VfreqMin
    END IF
  ELSEIF vfreq > VfreqMax THEN
     vfreq = VfreqMax
  END IF
  IF yres > ActiveLinesLimit AND Rotation = 0 THEN
    interlace = 2
    result = result OR %InterlacedRes
    IF yres < VirtualLinesLimit THEN ResVirtualize ( xres, yres, vfreq, hfreq ) : result = result OR %VirtualizedRes
  END IF
' here we start calculating an estimation of hfreq required for the mode, and check if it exceeds our monitor's limits.
  hfreq = vfreq * yres / ( interlace * ( 1 - vfreq * VerticalBlank ) )
  IF hfreq < HfreqMin THEN
    hfreq = HfreqMin
  ELSEIF hfreq > HfreqMax THEN
    IF yres > ActiveLinesLimit THEN
      ResVirtualize ( xres, yres, vfreq, hfreq )
      interlace = 2
      result = result OR %VirtualizedRes
    ELSE
      hfreq = HfreqMax
      VBlankLines = ROUND ( hfreq * VerticalBlank, 0 )
      vfreq = hfreq / ( yres / interlace + VBlankLines )
      result = result OR %BadVfreq
    END IF
  END IF
' Now, we need to know the minimum integer number of total lines needed to produce a mode that matches vfreq within
' our hfreq limits (vvtIni)
  vvtIni = ROUND ( hfreq / vfreq, 0 ) + IIF ( interlace = 2, 0.5, 0 )
  WHILE vfreq * vvtIni < HfreqMin - HfreqTolerance
    INCR vvtIni
  WEND
' Next block is used to get horizontal values and dotclock. This version uses an iterated loop and a dotclock table for
' accuracy, but it's not necessary and can be simplified like this (not tested):
'
'   hfreq = vfreq * vvtIni
'   ModelineGetLineParams ( hfreq, xres, Hhh, Hhi, Hhf, Hht )
'   DotClockReq = hfreq * Hht
'
' (notice hfreq is recalculated here from our integer vvtIni)
  FOR i = 0 TO Iterations
    hfreq = vfreq * ( vvtIni + i )
    IF hfreq <= HfreqMax + HfreqTolerance THEN
      ModelineGetLineParams ( hfreq, xres, newHhh, newHhi, newHhf, newHht )
      DotClockReq = hfreq * newHht
      ARRAY SCAN DotClockTable(), >= DotClockReq, TO j
      IF ABS ( DotClockTable( j - 1 ) - DotClockReq ) < ABS ( DotclockTable( j ) - DotClockReq ) THEN DECR j
      newVfreq = DotClockTable ( j ) / ( ( vvtIni + i ) * newHht )
      newDiff = ABS ( newVfreq - vfreq )
      IF newDiff < Diff OR Diff = 0 THEN
        hhh = newHhh : hhi = newHhi : hhf = NewHhf : hht = NewHht
        Diff = newDiff
        vIncr = i
        DotClockIdx = j
      END IF
    END IF
  NEXT
' Now we calculate the number of lines needed for vertical blanking. VerticalBlank is the time in µs our monitor needs
' to complete it's vertical retrace (front and back porches + sync pulses). I use 1280 µs for TVs and lowres monitors.
' Notice the 0.5 increment in case of interlace: this is used to make sure that vvt will be an uneven value if the mode
' is interlaced, otherwise our timings would be wrong as all interlaced modes have an uneven number of lines despite
' we set vvt as an even number.
  VBlankLines = INT ( hfreq * VerticalBlank ) + IIF ( interlace = 2, 0.5, 0 )
' Now we work out the vertical part of the modeline (forget about vIncr). We start with vvt and substract vvv and Vblanklines
' The rest is considered as fill in margins needed to achive the required vfreq, and we divide then by 2 (up and down margins
' are equal.
' For the sync pulse, 3 and 5 lines long pulses are hardcoded (progressive or interlaced) according to PAL standard, however
' this should be checked for higher frequencies
' Notice all calculations have beed done as 'progressive', and here we use the 'interlace' multiplier to get the right
' values in case of interlace
  vvv = yres
  vvt = ( vvtIni + vIncr ) * interlace
  margen = vvt - vvv - VBlankLines * interlace
  vvi = vvv + ROUND ( margen / 2, 0 ) + 1 * interlace
  vvf = vvi + IIF ( interlace = 2, 5, 3 )
' That's all, the rest is for storing values, etc.
  DotClock = DotClockIdx * 10
  hfreqReal = DotClockTable ( DotClockIdx ) / hht
  vfreqReal = hfreqReal / vvt * interlace
  vfreqLabel = ROUND ( vfreqReal, 0 )
  Mode(index).x = hhh
  Mode(index).y = vvv
  Mode(index).vfreq = vfreqReal
  Mode(index).Xlabel = hhh
  Mode(index).Ylabel = vvv
  Mode(index).Vlabel = vfreqLabel
  Mdln(index).dotclockReal = DotClockTable ( DotClockIdx )
  Mdln(index).dotclock = DotClock
  Mdln(index).hhh = hhh
  Mdln(index).hhi = hhi
  Mdln(index).hhf = hhf
  Mdln(index).hht = hht
  Mdln(index).vvv = vvv
  Mdln(index).vvi = vvi
  Mdln(index).vvf = vvf
  Mdln(index).vvt = vvt
  Mdln(index).interlace = interlace
  FUNCTION = result
END FUNCTION