SUBROUTINE cartfr(xlow,ylow,dx,dy,wl,hll,hbl, xc, yc, xm, ym, blen)
  IMPLICIT DOUBLE PRECISION (a-h,o-z)
  EXTERNAL fss
  LOGICAL fl(5),alll,allr
  DIMENSION x(10),y(10),xx(5),yy(5)
  COMMON/fsscorn/ xc0,yc0,xc1,yc1

!!
!!     # compute wl, fraction of cell that lies in left state.
!!     # For initial data with two states ql and qr separated by a
!!     # discontinuity. The curve along which the discontinuity lies is
!!     # specified by the function fdisc, which should return a value that
!!     # is negative on the side where ql lies and positive on the qr side.
!!
!!     # xlow,ylow is the coordinate of the lower left corner of the cell.
!!     # dx, dy are grid spacing in x and y.
!!
  xx(1) = xlow
  xx(2) = xlow
  xx(3) = xlow+dx
  xx(4) = xlow+dx
  xx(5) = xx(1)

  yy(1) = ylow
  yy(2) = ylow+dy
  yy(3) = ylow+dy
  yy(4) = ylow
  yy(5) = yy(1)
  alll = .TRUE.
  allr = .TRUE.

  zeps = 1e-12

  DO  i=1,4
     !! All in (-1) or all out (+1)
     fl(i) = fdisc(xx(i),yy(i)) .LT. 0.d0
     alll = alll .AND. fl(i)
     allr = allr .AND. (.NOT. fl(i))
  END DO

  fl(5) = fl(1)

  IF (alll) THEN
     wl = 1.d0
     hll = 1.d0
     hbl = 1.d0
     xm = xlow + dx/2
     ym = ylow + dy/2
     xc = xm
     yc = ym
     blen = 0
     RETURN
  ENDIF

  IF (allr) THEN
     wl = 0.d0
     hll = 0.d0
     hbl = 0.d0
     xm = xlow + dx/2
     ym = ylow + dy/2
     xc = xm
     yc = ym
     blen = 0
     RETURN
  ENDIF

  iv = 0
  DO  i=1,4
     IF (fl(i)) THEN
        iv = iv+1
        x(iv) = xx(i)
        y(iv) = yy(i)
     ENDIF
     IF (fl(i).NE.fl(i+1)) THEN
        iv = iv+1
        xc0 = xx(i)
        yc0 = yy(i)
        xc1 = xx(i+1)
        yc1 = yy(i+1)
        ss = zeroin(0.d0, 1.d0, fss, zeps)
        x(iv) = xx(i) + ss*(xx(i+1)-xx(i))
        y(iv) = yy(i) + ss*(yy(i+1)-yy(i))
     ENDIF

     IF (i.EQ.1) THEN
        IF (fl(i) .AND. fl(i+1)) THEN
           hll = 1.d0
           ym = y(i) + dy/2
        ELSE IF ((.NOT.fl(i)) .AND. (.NOT.fl(i+1))) THEN
           hll = 0.d0
           ym = 0
        ELSE IF (fl(i)) THEN
           hll = ss
           ym = y(i) + dy*ss/2
        ELSE
           hll = 1.d0 - ss
           ym = y(i) + dy*(1-ss)/2
        ENDIF
     ENDIF

     IF (i.EQ.4) THEN
        IF (fl(i) .AND. fl(i+1)) THEN
           hbl = 1.d0
           xm = x(i) + dx/2
        ELSE IF ((.NOT.fl(i)) .AND. (.NOT.fl(i+1))) THEN
           hbl = 0.d0
           xm = 0
        ELSE IF (fl(i)) THEN
           hbl = ss
           xm = x(i) - ss*dx/2  !! Assuming x(i) > x(i+1)
        ELSE
           hbl = 1.d0 - ss
           xm = x(i) - dx*(1-ss)/2
        ENDIF
     ENDIF
  END DO


!!     # compute area:

  x(iv+1) = x(1)
  y(iv+1) = y(1)
  area = 0.d0
  DO  i=1,iv
     area = area + .5d0*(y(i)+y(i+1))*(x(i+1)-x(i))
  END DO

  wl = area / (dx*dy)

  !! Compute center of mass for each polygon.
  IF (wl < 1e-12) THEN  !! Take average value if area too small
     DO i = 1,iv
        xc = xc + x(i)/iv
        yc = yc + y(i)/iv
     END DO
  ELSE
     sx = 0
     sy = 0
     DO i = 1,iv
        sx = sx - (y(i+1) - y(i))*(x(i+1)**2 + x(i+1)*x(i) + x(i)**2)/3.0
        sy = sy - (x(i+1) - x(i))*(y(i+1)**2 + y(i+1)*y(i) + y(i)**2)/3.0
     END DO
     xc = sx/(2.0*area)
     yc = -sy/(2.0*area)
  END IF

  !! Compute blen (length of short side)
  tol = 1e-15
  DO i = 1,iv
     xf = 0.5*(x(i+1) - x(i))/dx  !!
     yf = 0.5*(y(i+1) - y(i))/dy
     IF (ABS(xf) > tol .AND. ABS(xf) < (1-tol) .AND. &
          & (ABS(yf) > tol .AND. ABS(yf) < (1-tol))) THEN
        blen = SQRT((x(i+1) - x(i))**2 + (y(i+1) - y(i))**2)
        EXIT
     END IF
  END DO

  RETURN
END SUBROUTINE cartfr


DOUBLE PRECISION FUNCTION fss(s)
  IMPLICIT DOUBLE PRECISION (a-h,o-z)
  COMMON/fsscorn/ xc0,yc0,xc1,yc1

!! # compute fdisc at distance s between corners (xc0,yc0) and (xc1,yc1)

  x = xc0 + s*(xc1-xc0)
  y = yc0 + s*(yc1-yc0)
  fss = fdisc(x,y)
  RETURN
END FUNCTION fss

!! Special function required for cartfr.  Don't know why it was
!! written this way!
REAL FUNCTION fdisc(x,y)
  fdisc = -fbody(x,y)
END FUNCTION fdisc
