52:       double dsx = celldx[i];

53:       double dsy = celldy[j];

54:       double cc = soundspeed(i, j) * soundspeed(i, j);

55:       cc = cc + 2.0 * viscosity_a(i, j) / density0(i, j);

56:       cc = std::fmax(std::sqrt(cc), g_small);

57:       double dtct = dtc_safe * std::fmin(dsx, dsy) / cc;

58:       double div = 0.0;

59:       double dv1 = (xvel0(i, j) + xvel0(i + 0, j + 1)) * xarea(i, j);

60:       double dv2 = (xvel0(i + 1, j + 0) + xvel0(i + 1, j + 1)) * xarea(i + 1, j + 0);

61:       div = div + dv2 - dv1;

62:       double dtut = dtu_safe * 2.0 * volume(i, j) / std::fmax(std::fmax(std::fabs(dv1), std::fabs(dv2)), g_small * volume(i, j));

63:       dv1 = (yvel0(i, j) + yvel0(i + 1, j + 0)) * yarea(i, j);

64:       dv2 = (yvel0(i + 0, j + 1) + yvel0(i + 1, j + 1)) * yarea(i + 0, j + 1);

65:       div = div + dv2 - dv1;

66:       double dtvt = dtv_safe * 2.0 * volume(i, j) / std::fmax(std::fmax(std::fabs(dv1), std::fabs(dv2)), g_small * volume(i, j));

67:       div = div / (2.0 * volume(i, j));

68:       double dtdivt;

69:       if (div < -g_small) {

70:         dtdivt = dtdiv_safe * (-1.0 / div);

71:       } else {

72:         dtdivt = g_big;

73:       }

74:       double mins = std::fmin(dtct, std::fmin(dtut, std::fmin(dtvt, std::fmin(dtdivt, g_big))));

75:       dt_min_val0 = std::fmin(mins, dt_min_val0);