subroutine sw_compute_solve_diagnostics_adj(dt, nCells, nEdges, nVertices, &
     maxEdges, maxEdges2, nVertLevels, vertexDegree, &
     config_monotonic, config_thickness_adv_order, config_apvm_upwinding, &
     h, h_ad, u, u_ad, v_in, v_ad, vh, vh_ad, h_edge, h_edge_ad, &
     h_vertex_in, h_vertex_ad, circulation_in, circulation_ad, &
     vorticity_in, vorticity_ad, divergence, divergence_ad, ke, ke_ad, &
     pv_edge, pv_edge_ad, pv_vertex_in, pv_vertex_ad, pv_cell_in, pv_cell_ad, &
     vorticity_cell, vorticity_cell_ad,  gradPVn_in, gradPVn_ad, &
     gradPVt_in, gradPVt_ad, weightsOnEdge, kiteAreasOnVertex, cellsOnEdge, &
     cellsOnVertex, cellsOnCell, verticesOnEdge, nEdgesOnCell, edgesOnCell, &
     nEdgesOnEdge, edgesOnEdge, edgesOnVertex, dcEdge, dvEdge, areaCell, &
     areaTriangle, fVertex, deriv_two, boundaryEdge, boundaryCell)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  ! Compute diagnostic fields used in the tendency computations
  !
  ! Input: grid - grid metadata
  !
  ! Output: s - computed diagnostics
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

  ! ----- forward calculation requirement -----
  !needs updated values of h and u
  ! ----- end -----

  implicit none

  integer, parameter :: RKIND  = selected_real_kind(12)
  logical :: config_monotonic
  integer :: config_thickness_adv_order
  real (kind=RKIND) :: config_apvm_upwinding
  real (kind=RKIND), intent(in) :: dt
  integer :: nCells, nEdges, nVertices, nVertLevels, vertexDegree
  integer :: maxEdges, maxEdges2

  real(kind=RKIND), dimension(nVertices) :: fVertex
  real(kind=RKIND), dimension(nEdges) :: dvEdge
  real(kind=RKIND), dimension(nEdges) :: dcEdge
  real(kind=RKIND), dimension(nCells) :: areaCell
  real(kind=RKIND), dimension(nVertices) :: areaTriangle
  real(kind=RKIND), dimension(nVertLevels,nEdges) :: vh
  real(kind=RKIND), dimension(nVertLevels,nEdges) :: vh_ad
  real(kind=RKIND), dimension(maxEdges2,nEdges) :: weightsOnEdge
  real(kind=RKIND), dimension(vertexDegree,nVertices) :: kiteAreasOnVertex
  real(kind=RKIND), dimension(nVertLevels,nEdges) :: h_edge
  real(kind=RKIND), dimension(nVertLevels,nEdges) :: h_edge_ad
  real(kind=RKIND), dimension(nVertLevels,nCells), intent(in) :: h
  real(kind=RKIND), dimension(nVertLevels,nCells) :: h_ad
  real(kind=RKIND), dimension(nVertLevels,nEdges), intent(in) :: u
  real(kind=RKIND), dimension(nVertLevels,nEdges) :: u_ad
  real(kind=RKIND), dimension(nVertLevels,nEdges) :: v_in
  real(kind=RKIND), dimension(nVertLevels,nEdges) :: v
  real(kind=RKIND), dimension(nVertLevels,nEdges) :: v_ad
  real(kind=RKIND), dimension(nVertLevels,nVertices) :: circulation_in
  real(kind=RKIND), dimension(nVertLevels,nVertices) :: circulation
  real(kind=RKIND), dimension(nVertLevels,nVertices) :: circulation_ad
  real(kind=RKIND), dimension(nVertLevels,nVertices) :: vorticity_in
  real(kind=RKIND), dimension(nVertLevels,nVertices) :: vorticity
  real(kind=RKIND), dimension(nVertLevels,nVertices) :: vorticity_ad
  real(kind=RKIND), dimension(nVertLevels,nCells) :: ke
  real(kind=RKIND), dimension(nVertLevels,nCells) :: ke_ad
  real(kind=RKIND), dimension(nVertLevels,nEdges) :: pv_edge
  real(kind=RKIND), dimension(nVertLevels,nEdges) :: pv_edge_ad
  real(kind=RKIND), dimension(nVertLevels,nVertices) :: pv_vertex_in
  real(kind=RKIND), dimension(nVertLevels,nVertices) :: pv_vertex
  real(kind=RKIND), dimension(nVertLevels,nVertices) :: pv_vertex_ad
  real(kind=RKIND), dimension(nVertLevels,nCells) :: pv_cell_in
  real(kind=RKIND), dimension(nVertLevels,nCells) :: pv_cell
  real(kind=RKIND), dimension(nVertLevels,nCells) :: pv_cell_ad
  real(kind=RKIND), dimension(nVertLevels,nEdges) :: gradPVn_in
  real(kind=RKIND), dimension(nVertLevels,nEdges) :: gradPVn
  real(kind=RKIND), dimension(nVertLevels,nEdges) :: gradPVn_ad
  real(kind=RKIND), dimension(nVertLevels,nEdges) :: gradPVt_in
  real(kind=RKIND), dimension(nVertLevels,nEdges) :: gradPVt
  real(kind=RKIND), dimension(nVertLevels,nEdges) :: gradPVt_ad
  real(kind=RKIND), dimension(nVertLevels,nCells) :: divergence
  real(kind=RKIND), dimension(nVertLevels,nCells) :: divergence_ad
  real(kind=RKIND), dimension(nVertLevels,nVertices) :: h_vertex_in
  real(kind=RKIND), dimension(nVertLevels,nVertices) :: h_vertex
  real(kind=RKIND), dimension(nVertLevels,nVertices) :: h_vertex_ad
  real(kind=RKIND), dimension(nVertLevels,nCells) :: vorticity_cell
  real(kind=RKIND), dimension(nVertLevels,nCells) :: vorticity_cell_ad

  integer, dimension(2,nEdges) :: cellsOnEdge
  integer, dimension(vertexDegree,nVertices) :: cellsOnVertex
  integer, dimension(2,nEdges) :: verticesOnEdge
  integer, dimension(maxEdges,nCells) :: edgesOnCell
  integer, dimension(maxEdges2,nEdges) :: edgesOnEdge
  integer, dimension(vertexDegree,nVertices) :: edgesOnVertex
  integer, dimension(nVertLevels,nEdges) :: boundaryEdge
  integer, dimension(nVertLevels,nCells) :: boundaryCell
  integer, dimension(maxEdges,nCells) :: cellsOnCell

  integer, dimension(nCells) :: nEdgesOnCell
  integer, dimension(nEdges) :: nEdgesOnEdge
  real (kind=RKIND), dimension(15,2,nEdges) :: deriv_two

  ! ----- local variables -----
  real (kind=RKIND) :: r, h1, h2
  real (kind=RKIND) :: d2fdx2_cell1, d2fdx2_cell2
  real (kind=RKIND) :: d2fdx2_cell1_ad, d2fdx2_cell2_ad
  real (kind=RKIND) :: coef_3rd_order

  integer :: iEdge, iCell, iVertex, k, cell1, cell2, vertex1, vertex2, eoe, i, j, cov
  real (kind=RKIND) :: flux, vorticity_abs, workpv
  real (kind=RKIND) :: flux_ad, vorticity_abs_ad, workpv_ad
  ! ----- end local variables -----

  ! Find those cells that have an edge on the boundary
  boundaryCell(:,:) = 0
  do iEdge = 1, nEdges
     do k = 1, nVertLevels
        if(boundaryEdge(k,iEdge)==1) then
           cell1 = cellsOnEdge(1,iEdge)
           cell2 = cellsOnEdge(2,iEdge)
           boundaryCell(k,cell1) = 1
           boundaryCell(k,cell2) = 1
        endif
     enddo
  enddo

  coef_3rd_order = 0.
  if (config_thickness_adv_order == 3) coef_3rd_order = 1.0
  if (config_thickness_adv_order == 3 .and. config_monotonic) coef_3rd_order = 0.25

  ! ----- forward calculation -----
  circulation(:,:) = 0.0
  do iEdge = 1, nEdges
     do k = 1, nVertLevels
        circulation(k,verticesOnEdge(1,iEdge)) = circulation(k,verticesOnEdge(1,iEdge)) - dcEdge(iEdge) * u(k,iEdge)
        circulation(k,verticesOnEdge(2,iEdge)) = circulation(k,verticesOnEdge(2,iEdge)) + dcEdge(iEdge) * u(k,iEdge)
     end do
  end do
  do iVertex = 1, nVertices
     do k = 1, nVertLevels
        vorticity(k,iVertex) = circulation(k,iVertex) / areaTriangle(iVertex)
     end do
  end do

  v(:,:) = 0.0
  do iEdge = 1,nEdges
     do i = 1, nEdgesOnEdge(iEdge)
        eoe = edgesOnEdge(i,iEdge)
        do k = 1,nVertLevels
           v(k,iEdge) = v(k,iEdge) + weightsOnEdge(i,iEdge) * u(k, eoe)
        end do
     end do
  end do

  do iVertex = 1,nVertices
     do k = 1, nVertLevels
        h_vertex(k,iVertex) = 0.0
        do i = 1, vertexDegree
           h_vertex(k,iVertex) = h_vertex(k,iVertex) + h(k,cellsOnVertex(i,iVertex)) * kiteAreasOnVertex(i,iVertex)
        end do
        h_vertex(k,iVertex) = h_vertex(k,iVertex) / areaTriangle(iVertex)
        pv_vertex(k,iVertex) = (fVertex(iVertex) + vorticity(k,iVertex)) / h_vertex(k,iVertex)
     end do
  end do

  do iEdge = 1, nEdges
     do k = 1, nVertLevels
        gradPVt(k,iEdge) = (pv_vertex(k,verticesOnEdge(2,iEdge)) - pv_vertex(k,verticesOnEdge(1,iEdge))) / dvEdge(iEdge)
     enddo
  enddo

  pv_cell(:,:) = 0.0
  do iVertex = 1, nVertices
     do i = 1, vertexDegree
        iCell = cellsOnVertex(i,iVertex)
        if (iCell <= nCells) then
           do k = 1, nVertLevels
              pv_cell(k,iCell) = pv_cell(k,iCell) &
                   + kiteAreasOnVertex(i, iVertex) * pv_vertex(k, iVertex) &
                   / areaCell(iCell)
           end do
        end if
     end do
  end do

  gradPVn(:,:) = 0.0
  do iEdge = 1, nEdges
     if( cellsOnEdge(1,iEdge) <= nCells .and. cellsOnEdge(2,iEdge) <= nCells) then
        do k = 1, nVertLevels
           gradPVn(k,iEdge) = (pv_cell(k,cellsOnEdge(2,iEdge)) - pv_cell(k,cellsOnEdge(1,iEdge))) / dcEdge(iEdge)
        enddo
     endif
  enddo

  ! ----- end forward calculation -----

  ! Modify PV edge with upstream bias.
  do iEdge = nEdges, 1, -1
     do k = nVertLevels, 1, -1
        gradPVn_ad(k,iEdge)=gradPVn_ad(k,iEdge) &
             -config_apvm_upwinding * u(k,iEdge) * dt * pv_edge_ad(k,iEdge)
        u_ad(k,iEdge)=u_ad(k,iEdge) &
             -config_apvm_upwinding * dt * gradPVn(k,iEdge) * pv_edge_ad(k,iEdge)
     enddo
  enddo

  do iEdge = nEdges, 1, -1
     if( cellsOnEdge(1,iEdge) <= nCells .and. cellsOnEdge(2,iEdge) <= nCells) then
        do k = nVertLevels, 1, -1
           pv_cell_ad(k,cellsOnEdge(1,iEdge))=pv_cell_ad(k,cellsOnEdge(1,iEdge)) &
                -gradPVn_ad(k,iEdge)/dcEdge(iEdge)
           pv_cell_ad(k,cellsOnEdge(2,iEdge))=pv_cell_ad(k,cellsOnEdge(2,iEdge)) &
                +gradPVn_ad(k,iEdge)/dcEdge(iEdge)
           gradPVn_ad(k,iEdge)=0.
        enddo
     endif
  enddo
  gradPVn_ad(:,:) = 0.0

  do iVertex = nVertices, 1, -1
     do i = vertexDegree, 1, -1
        iCell = cellsOnVertex(i,iVertex)
        if (iCell <= nCells) then
           do k = nVertLevels, 1, -1
              vorticity_ad(k, iVertex)=vorticity_ad(k, iVertex) &
                   +kiteAreasOnVertex(i, iVertex)/areaCell(iCell) * vorticity_cell_ad(k,iCell)
              pv_vertex_ad(k, iVertex)=pv_vertex_ad(k, iVertex) &
                   +kiteAreasOnVertex(i, iVertex)/areaCell(iCell) * pv_cell_ad(k,iCell)
           enddo
        endif
     enddo
  enddo
  vorticity_cell_ad(:,:) = 0.0
  pv_cell_ad(:,:) = 0.0

  ! ----- forward calculation -----
  !gradPVt=???
  !v=???
  ! ----- end forward calculation -----
  do iEdge = nEdges, 1, -1
     do k = nVertLevels, 1, -1
        gradPVt_ad(k,iEdge)=gradPVt_ad(k,iEdge) &
             -config_apvm_upwinding * v(k,iEdge) * dt * pv_edge_ad(k,iEdge)
        v_ad(k,iEdge)=v_ad(k,iEdge) &
             -config_apvm_upwinding * dt * gradPVt(k,iEdge) * pv_edge_ad(k,iEdge)
     end do
  end do

  do iVertex = nVertices, 1, -1
     do i = vertexDegree, 1, -1
        iEdge = edgesOnVertex(i,iVertex)
        do k = nVertLevels, 1, -1
           pv_vertex_ad(k,iVertex)=pv_vertex_ad(k,iVertex) + 0.5*pv_edge_ad(k,iEdge)
        end do
     end do
  end do
  pv_edge_ad(:,:) = 0.0

  do iEdge = nEdges, 1, -1
     do k = nVertLevels, 1, -1
        pv_vertex_ad(k,verticesOnEdge(1,iEdge))=pv_vertex_ad(k,verticesOnEdge(1,iEdge)) &
             -gradPVt_ad(k,iEdge)/dvEdge(iEdge)
        pv_vertex_ad(k,verticesOnEdge(2,iEdge))=pv_vertex_ad(k,verticesOnEdge(2,iEdge)) &
             +gradPVt_ad(k,iEdge)/dvEdge(iEdge)
        gradPVt_ad(k,iEdge)=0.
     end do
  end do

  ! ----- forward calculation -----
  !vorticity=???
  !h_vertex=???
  ! ----- end forward calculation -----
  do iVertex = nVertices, 1, -1
     do k = nVertLevels, 1, -1
        ! ----- forward calculation -----
        !h_vertex(k,iVertex) = 0.0
        !do i = 1, vertexDegree
        !   h_vertex(k,iVertex) = h_vertex(k,iVertex) + h(k,cellsOnVertex(i,iVertex)) * kiteAreasOnVertex(i,iVertex)
        !end do
        !h_vertex(k,iVertex) = h_vertex(k,iVertex) / areaTriangle(iVertex)
        ! ----- end forward calculation -----
        h_vertex_ad(k,iVertex)=h_vertex_ad(k,iVertex) &
             -(fVertex(iVertex) + vorticity(k,iVertex))/(h_vertex(k,iVertex)**2)&
             *pv_vertex_ad(k,iVertex)
        vorticity_ad(k,iVertex)=vorticity_ad(k,iVertex) &
             +pv_vertex_ad(k,iVertex)/h_vertex(k,iVertex)
        pv_vertex_ad(k,iVertex)=0.

        h_vertex_ad(k,iVertex)=h_vertex_ad(k,iVertex)/areaTriangle(iVertex)

        do i = vertexDegree, 1, -1
           h_ad(k,cellsOnVertex(i,iVertex))=h_ad(k,cellsOnVertex(i,iVertex)) &
                +kiteAreasOnVertex(i,iVertex)*h_vertex_ad(k,iVertex)
        end do
        h_vertex_ad(k,iVertex) = 0.0
     end do
  end do

  ! vh omitted

  do iEdge = nEdges, 1, -1
     do i = nEdgesOnEdge(iEdge), 1, -1
        eoe = edgesOnEdge(i,iEdge)
        do k = nVertLevels, 1, -1
           u_ad(k, eoe)=u_ad(k, eoe) + weightsOnEdge(i,iEdge)*v_ad(k,iEdge)
        end do
     end do
  end do
  v_ad(:,:)=0.

  do iCell = nCells, 1, -1
     do k=nVertLevels, 1, -1
        ke_ad(k,iCell)=ke_ad(k,iCell)/areaCell(iCell)
     end do
     do i = nEdgesOnCell(iCell), 1, -1
        iEdge = edgesOnCell(i,iCell)
        do k = nVertLevels, 1, -1
           u_ad(k,iEdge)=u_ad(k,iEdge) &
                +0.25*dcEdge(iEdge)*dvEdge(iEdge)*2.*u(k,iEdge)*ke_ad(k,iCell)
        end do
     end do
  end do
  ke_ad(:,:)=0.

  do iCell = nCells, 1, -1
     r = 1.0 / areaCell(iCell)
     do k = nVertLevels, 1, -1
        divergence_ad(k,iCell)=divergence_ad(k,iCell)*r
     end do
  end do
  do iEdge = nEdges, 1, -1
     cell1 = cellsOnEdge(1,iEdge)
     cell2 = cellsOnEdge(2,iEdge)
     if(cell2 <= nCells) then
        do k = nVertLevels, 1, -1
           u_ad(k,iEdge)=u_ad(k,iEdge) - dvEdge(iEdge)*divergence_ad(k,cell2)
        end do
     end if
     if(cell1 <= nCells) then
        do k = nVertLevels, 1, -1
           u_ad(k,iEdge)=u_ad(k,iEdge) + dvEdge(iEdge)*divergence_ad(k,cell1)
        end do
     end if
  end do
  divergence_ad(:,:) = 0.0

  do iVertex = nVertices, 1, -1
     do k = nVertLevels, 1, -1
        circulation_ad(k,iVertex)=circulation_ad(k,iVertex) &
          + vorticity_ad(k,iVertex)/areaTriangle(iVertex)
        vorticity_ad(k,iVertex)=0.
     end do
  end do
  do iEdge = nEdges, 1, -1
     do k = nVertLevels, 1, -1
        u_ad(k,iEdge)=u_ad(k,iEdge) + dcEdge(iEdge)*circulation_ad(k,verticesOnEdge(2,iEdge))
        u_ad(k,iEdge)=u_ad(k,iEdge) - dcEdge(iEdge)*circulation_ad(k,verticesOnEdge(1,iEdge))
     end do
  end do
  circulation_ad(:,:) = 0.0

  if (config_thickness_adv_order == 2) then
     do iEdge = nEdges, 1, -1
        cell1 = cellsOnEdge(1,iEdge)
        cell2 = cellsOnEdge(2,iEdge)
        if (cell1 <= nCells .and. cell2 <= nCells) then
           do k = nVertLevels, 1, -1
              h_ad(k,cell2)=h_ad(k,cell2) + 0.5*h_edge_ad(k,iEdge)
              h_ad(k,cell1)=h_ad(k,cell1) + 0.5*h_edge_ad(k,iEdge)
              h_edge_ad(k,iEdge)=0.
           end do
        end if
     end do
  else if (config_thickness_adv_order == 3) then
     do iEdge = nEdges, 1, -1
        cell1 = cellsOnEdge(1,iEdge)
        cell2 = cellsOnEdge(2,iEdge)
        if (cell1 <= nCells .and. cell2 <= nCells) then
           do k = 1, nVertLevels
              if (u(k,iEdge) > 0) then
                 d2fdx2_cell2_ad= (dcEdge(iEdge)**2)*coef_3rd_order/12.*h_edge_ad(k,iEdge)
                 d2fdx2_cell1_ad=-(dcEdge(iEdge)**2)*coef_3rd_order/12.*h_edge_ad(k,iEdge)
                 d2fdx2_cell2_ad= d2fdx2_cell2_ad &
                      -(dcEdge(iEdge)**2)/12.*h_edge_ad(k,iEdge)
                 d2fdx2_cell1_ad= d2fdx2_cell1_ad &
                      -(dcEdge(iEdge)**2)/12.*h_edge_ad(k,iEdge)
                 h_ad(k,cell2)=h_ad(k,cell2) + 0.5*h_edge_ad(k,iEdge)
                 h_ad(k,cell1)=h_ad(k,cell1) + 0.5*h_edge_ad(k,iEdge)
                 h_edge_ad(k,iEdge)=0.
              else
                 d2fdx2_cell2_ad=-(dcEdge(iEdge)**2)*coef_3rd_order/12.*h_edge_ad(k,iEdge)
                 d2fdx2_cell1_ad= (dcEdge(iEdge)**2)*coef_3rd_order/12.*h_edge_ad(k,iEdge)
                 d2fdx2_cell2_ad=d2fdx2_cell2_ad &
                      -(dcEdge(iEdge)**2)/12.*h_edge_ad(k,iEdge)
                 d2fdx2_cell1_ad=d2fdx2_cell1_ad &
                      -(dcEdge(iEdge)**2)/12.*h_edge_ad(k,iEdge)
                 h_ad(k,cell2)=h_ad(k,cell2) + 0.5*h_edge_ad(k,iEdge)
                 h_ad(k,cell1)=h_ad(k,cell1) + 0.5*h_edge_ad(k,iEdge)
                 h_edge_ad(k,iEdge)=0.
              end if

              if(boundaryCell(k,cell1) == 0 .and. boundaryCell(k,cell2) == 0) then
                 do i=nEdgesOnCell(cell2), 1, -1
                    h_ad(k, cellsOnCell(i,cell2))=h_ad(k, cellsOnCell(i,cell2)) &
                         +deriv_two(i+1,2,iEdge)*d2fdx2_cell2_ad
                 end do
                 do i = nEdgesOnCell(cell1), 1, -1
                    h_ad(k, cellsOnCell(i,cell1))=h_ad(k, cellsOnCell(i,cell1)) &
                         +deriv_two(i+1,1,iEdge)*d2fdx2_cell1_ad
                 end do
                 h_ad(k,cell2)=h_ad(k,cell2) + deriv_two(1,2,iEdge)*d2fdx2_cell2_ad
                 d2fdx2_cell2_ad=0.
                 h_ad(k,cell1)=h_ad(k,cell1) + deriv_two(1,1,iEdge)*d2fdx2_cell1_ad
                 d2fdx2_cell1_ad=0.
              end if ! (boundaryCell(k,cell1) .eq. 0 .and.
              d2fdx2_cell2_ad = 0.0
              d2fdx2_cell1_ad = 0.0
           end do ! k = 1, nVertLevels
        end if ! (cell1 <= nCells .and. cell2 <= nCells) then
     end do ! iEdge = nEdges, 1, -1

  else if (config_thickness_adv_order == 4) then
     do iEdge = nEdges, 1, -1
        cell1 = cellsOnEdge(1,iEdge)
        cell2 = cellsOnEdge(2,iEdge)

        if (cell1 <= nCells .and. cell2 <= nCells) then
           do k = nVertLevels, 1, -1
              d2fdx2_cell2_ad=-(dcEdge(iEdge) **2)/12. * h_edge_ad(k,iEdge)
              d2fdx2_cell1_ad=-(dcEdge(iEdge) **2)/12. * h_edge_ad(k,iEdge)
              h_ad(k,cell2)=h_ad(k,cell2) + 0.5*h_edge_ad(k,iEdge)
              h_ad(k,cell1)=h_ad(k,cell1) + 0.5*h_edge_ad(k,iEdge)
              h_edge_ad(k,iEdge)=0.

              if (boundaryCell(k,cell1) == 0 .and. boundaryCell(k,cell2) == 0) then
                 do i = nEdgesOnCell(cell2), 1, -1
                    h_ad(k, cellsOnCell(i,cell2))=h_ad(k, cellsOnCell(i,cell2)) &
                         +deriv_two(i+1,2,iEdge)*d2fdx2_cell2_ad
                 end do
                 do i = nEdgesOnCell(cell1), 1, -1
                    h_ad(k, cellsOnCell(i,cell1))=h_ad(k, cellsOnCell(i,cell1)) &
                         +deriv_two(i+1,1,iEdge)*d2fdx2_cell1_ad
                 end do
                 h_ad(k,cell2)=h_ad(k,cell2) + deriv_two(1,2,iEdge)*d2fdx2_cell2_ad
                 d2fdx2_cell2_ad=0.
                 h_ad(k,cell1)=h_ad(k,cell1) + deriv_two(1,1,iEdge)*d2fdx2_cell1_ad
                 d2fdx2_cell1_ad=0.
              end if ! (boundaryCell(k,cell1) .eq. 0 .and.
              d2fdx2_cell2_ad = 0.0
              d2fdx2_cell1_ad = 0.0
           end do ! k=nVertLevels, 1, -1
        end if ! (cell1 <= nCells .and. cell2 <= nCells) then
     end do ! iEdge = 1, nEdges

  end if ! (config_thickness_adv_order == 2) then

end subroutine sw_compute_solve_diagnostics_adj

subroutine sw_compute_tend_adj(nCells, nEdges, nVertices, nVertLevels, &
     maxEdges, maxEdges2, vertexDegree, &
     config_bottom_drag, config_wind_stress, config_h_mom_eddy_visc2, &
     config_h_mom_eddy_visc4, h, h_ad, u, u_ad, v, v_ad, h_edge, h_edge_ad, &
     circulation, circulation_ad, vorticity, vorticity_ad, &
     divergence, divergence_ad, ke, ke_ad, pv_edge, pv_edge_ad, vh, vh_ad, &
     weightsOnEdge, kiteAreasOnVertex, cellsOnEdge, cellsOnVertex, &
     verticesOnEdge, nEdgesOnCell, edgesOnCell, nEdgesOnEdge, &
     edgesOnEdge, edgesOnVertex, dcEdge, dvEdge, areaCell, areaTriangle, &
     h_s, fVertex, fEdge, u_src, meshScalingDel2, meshScalingDel4, &
     tend_h, tend_h_ad, tend_u, tend_u_ad)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  ! Compute height and normal wind tendencies, as well as diagnostic variables
  ! Input: s - current model state
  !        grid - grid metadata
  ! Output: tend - computed tendencies for prognostic variables
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

  ! ----- forward calculation requirement -----
  ! needs updated values of ke, h_edge, and pv_edge
  ! ----- end -----
  implicit none

  integer, parameter :: RKIND  = selected_real_kind(12)
  real (kind=RKIND), parameter :: gravity = 9.80616
  integer :: nCells, nEdges, nVertices, nVertLevels
  integer :: maxEdges, maxEdges2, vertexDegree
  real(kind=RKIND), dimension(nCells) :: h_s
  real(kind=RKIND), dimension(nVertices) :: fVertex
  real(kind=RKIND), dimension(nEdges) :: fEdge
  real(kind=RKIND), dimension(nEdges) :: dvEdge
  real(kind=RKIND), dimension(nEdges) :: dcEdge
  real(kind=RKIND), dimension(nCells) :: areaCell
  real(kind=RKIND), dimension(nVertices) :: areaTriangle
  real(kind=RKIND), dimension(nEdges) :: meshScalingDel2
  real(kind=RKIND), dimension(nEdges) :: meshScalingDel4
  real(kind=RKIND), dimension(nVertLevels,nEdges) :: vh
  real(kind=RKIND), dimension(nVertLevels,nEdges) :: vh_ad
  real(kind=RKIND), dimension(maxEdges2,nEdges) :: weightsOnEdge
  real(kind=RKIND), dimension(vertexDegree,nVertices) :: kiteAreasOnVertex
  real(kind=RKIND), dimension(nVertLevels,nEdges) :: h_edge
  real(kind=RKIND), dimension(nVertLevels,nEdges) :: h_edge_ad
  real(kind=RKIND), dimension(nVertLevels,nCells) :: h
  real(kind=RKIND), dimension(nVertLevels,nCells) :: h_ad
  real(kind=RKIND), dimension(nVertLevels,nEdges) :: u
  real(kind=RKIND), dimension(nVertLevels,nEdges) :: u_ad
  real(kind=RKIND), dimension(nVertLevels,nEdges) :: v
  real(kind=RKIND), dimension(nVertLevels,nEdges) :: v_ad
  real(kind=RKIND), dimension(nVertLevels,nCells) :: tend_h
  real(kind=RKIND), dimension(nVertLevels,nCells) :: tend_h_ad
  real(kind=RKIND), dimension(nVertLevels,nEdges) :: tend_u
  real(kind=RKIND), dimension(nVertLevels,nEdges) :: tend_u_ad
  real(kind=RKIND), dimension(nVertLevels,nVertices) :: circulation
  real(kind=RKIND), dimension(nVertLevels,nVertices) :: circulation_ad
  real(kind=RKIND), dimension(nVertLevels,nVertices) :: vorticity
  real(kind=RKIND), dimension(nVertLevels,nVertices) :: vorticity_ad
  real(kind=RKIND), dimension(nVertLevels,nCells) :: ke
  real(kind=RKIND), dimension(nVertLevels,nCells) :: ke_ad
  real(kind=RKIND), dimension(nVertLevels,nEdges) :: pv_edge
  real(kind=RKIND), dimension(nVertLevels,nEdges) :: pv_edge_ad
  real(kind=RKIND), dimension(nVertLevels,nCells) :: divergence
  real(kind=RKIND), dimension(nVertLevels,nCells) :: divergence_ad
  real(kind=RKIND), dimension(nVertLevels,nVertices) :: h_vertex
  real(kind=RKIND), dimension(nVertLevels,nVertices) :: h_vertex_ad

  integer, dimension(2,nEdges) :: cellsOnEdge
  integer, dimension(vertexDegree,nVertices) :: cellsOnVertex
  integer, dimension(2,nEdges) :: verticesOnEdge
  integer, dimension(maxEdges,nCells) :: edgesOnCell
  integer, dimension(maxEdges2,nEdges) :: edgesOnEdge
  integer, dimension(vertexDegree,nVertices) :: edgesOnVertex
  integer, dimension(nCells) :: nEdgesOnCell
  integer, dimension(nEdges) :: nEdgesOnEdge

  real (kind=RKIND), dimension(nVertLevels, nEdges) :: u_src

  logical :: config_wind_stress, config_bottom_drag
  real (kind=RKIND) :: config_h_mom_eddy_visc2, config_h_mom_eddy_visc4

  ! ----- local variables -----
  integer :: iEdge, iCell, iVertex, k, cell1, cell2, vertex1, vertex2, eoe, i, j, timeLevel
  real (kind=RKIND) :: flux, vorticity_abs, workpv, q, upstream_bias
  real (kind=RKIND) :: flux_ad, vorticity_abs_ad, workpv_ad, q_ad, upstream_bias_ad
  real (kind=RKIND) :: ke_edge, ke_edge_ad
  real (kind=RKIND), parameter :: rho_ref = 1000.0
  real (kind=RKIND) :: r, u_diffusion, u_diffusion_ad
  real (kind=RKIND), allocatable, dimension(:,:) :: delsq_divergence
  real (kind=RKIND), allocatable, dimension(:,:) :: delsq_divergence_ad
  real (kind=RKIND), allocatable, dimension(:,:) :: delsq_u
  real (kind=RKIND), allocatable, dimension(:,:) :: delsq_u_ad
  real (kind=RKIND), allocatable, dimension(:,:) :: delsq_circulation, delsq_vorticity
  real (kind=RKIND), allocatable, dimension(:,:) :: delsq_circulation_ad, delsq_vorticity_ad
  ! ----- end local variables

  if (config_bottom_drag) then
     do iEdge = nEdges, 1, -1
        ! ----- forward calculation -----
        ke_edge = 0.5 * ( ke(1,cellsOnEdge(1,iEdge)) + ke(1,cellsOnEdge(2,iEdge)))
        ! ----- end forward calculation -----

        h_edge_ad(1,iEdge)=h_edge_ad(1,iEdge) &
             +1.0e-3*u(1,iEdge)*sqrt(2.0*ke_edge)/(h_edge(1,iEdge)**2)*tend_u_ad(1,iEdge)
        ke_edge_ad=-0.5*1.0e-3*u(1,iEdge)*2.0/sqrt(2.0*ke_edge)/h_edge(1,iEdge)*tend_u_ad(1,iEdge)
        u_ad(1,iEdge)=u_ad(1,iEdge) &
             -1.0e-3*sqrt(2.0*ke_edge)/h_edge(1,iEdge)*tend_u_ad(1,iEdge)

        ke_ad(1,cellsOnEdge(2,iEdge))=ke_ad(1,cellsOnEdge(2,iEdge)) +0.5*ke_edge_ad
        ke_ad(1,cellsOnEdge(1,iEdge))=ke_ad(1,cellsOnEdge(1,iEdge)) +0.5*ke_edge_ad
        ke_edge_ad=0.
     end do
  end if

  if(config_wind_stress) then
     do iEdge = nEdges, 1, -1
        h_edge_ad(1,iEdge)=h_edge_ad(1,iEdge) &
             -u_src(1,iEdge) / rho_ref / (h_edge(1,iEdge)**2) * tend_u_ad(1,iEdge)
     end do
  end if

  if (config_h_mom_eddy_visc4 > 0.0) then
     allocate(delsq_divergence_ad(nVertLevels, nCells+1))
     allocate(delsq_divergence(nVertLevels, nCells+1))
     allocate(delsq_u_ad(nVertLevels, nEdges+1))
     allocate(delsq_u(nVertLevels, nEdges+1))
     allocate(delsq_circulation_ad(nVertLevels, nVertices+1))
     allocate(delsq_circulation(nVertLevels, nVertices+1))
     allocate(delsq_vorticity_ad(nVertLevels, nVertices+1))
     allocate(delsq_vorticity(nVertLevels, nVertices+1))
     delsq_divergence_ad =0.
     delsq_u_ad          =0.
     delsq_circulation_ad=0.
     delsq_vorticity_ad  =0.

     do iEdge = nEdges, 1, -1
        cell1 = cellsOnEdge(1,iEdge)
        cell2 = cellsOnEdge(2,iEdge)
        vertex1 = verticesOnEdge(1,iEdge)
        vertex2 = verticesOnEdge(2,iEdge)

        do k=nVertLevels, 1, -1
           u_diffusion_ad=-tend_u_ad(k,iEdge)
           u_diffusion_ad=meshScalingDel4(iEdge)*config_h_mom_eddy_visc4 * u_diffusion_ad

           delsq_vorticity_ad(k,vertex1)=delsq_vorticity_ad(k,vertex1) &
                +u_diffusion_ad/dvEdge(iEdge)
           delsq_vorticity_ad(k,vertex2)=delsq_vorticity_ad(k,vertex2) &
                -u_diffusion_ad/dvEdge(iEdge)
           delsq_divergence_ad(k,cell1)=delsq_divergence_ad(k,cell1) &
                -u_diffusion_ad/dcEdge(iEdge)
           delsq_divergence_ad(k,cell2)=delsq_divergence_ad(k,cell2) &
                +u_diffusion_ad/dcEdge(iEdge)
           u_diffusion_ad=0.
        end do
     end do ! iEdge = nEdges, 1, -1

     do iCell = nCells, 1, -1
        r = 1.0 / areaCell(iCell)
        do k = nVertLevels, 1, -1
           delsq_divergence_ad(k,iCell)=delsq_divergence_ad(k,iCell)*r
        end do
     end do

     do iEdge = nEdges, 1, -1
        cell1 = cellsOnEdge(1,iEdge)
        cell2 = cellsOnEdge(2,iEdge)
        do k=nVertLevels, 1, -1
           delsq_u_ad(k,iEdge)=delsq_u_ad(k,iEdge) -dvEdge(iEdge)*delsq_divergence_ad(k,cell2)
           delsq_u_ad(k,iEdge)=delsq_u_ad(k,iEdge) +dvEdge(iEdge)*delsq_divergence_ad(k,cell1)
        end do
     end do
     delsq_divergence_ad(:,:) = 0.0

     do iVertex = nVertices, 1, -1
        r = 1.0 / areaTriangle(iVertex)
        do k = nVertLevels, 1, -1
           delsq_circulation_ad(k,iVertex)=delsq_circulation_ad(k,iVertex) &
                +delsq_vorticity_ad(k,iVertex) * r
           delsq_vorticity_ad(k,iVertex)=0.
        end do
     end do

     do iEdge = nEdges, 1, -1
        vertex1 = verticesOnEdge(1,iEdge)
        vertex2 = verticesOnEdge(2,iEdge)
        do k=nVertLevels, 1, -1
           delsq_u_ad(k,iEdge)=delsq_u_ad(k,iEdge) + dcEdge(iEdge)*delsq_circulation_ad(k,vertex2)
           delsq_u_ad(k,iEdge)=delsq_u_ad(k,iEdge) - dcEdge(iEdge)*delsq_circulation_ad(k,vertex1)
        end do
     end do
     delsq_circulation_ad(:,:) = 0.0

     do iEdge = nEdges, 1, -1
        cell1 = cellsOnEdge(1,iEdge)
        cell2 = cellsOnEdge(2,iEdge)
        vertex1 = verticesOnEdge(1,iEdge)
        vertex2 = verticesOnEdge(2,iEdge)
        do k = nVertLevels, 1, -1
           vorticity_ad(k,vertex1)=vorticity_ad(k,vertex1)+delsq_u_ad(k,iEdge)/dvEdge(iEdge)
           vorticity_ad(k,vertex2)=vorticity_ad(k,vertex2)-delsq_u_ad(k,iEdge)/dvEdge(iEdge)

           divergence_ad(k,cell1)=divergence_ad(k,cell1)-delsq_u_ad(k,iEdge)/dcEdge(iEdge)
           divergence_ad(k,cell2)=divergence_ad(k,cell2)+delsq_u_ad(k,iEdge)/dcEdge(iEdge)
        end do
     end do
     delsq_u_ad(:,:) = 0.0

     deallocate(delsq_divergence_ad)
     deallocate(delsq_divergence)
     deallocate(delsq_u_ad)
     deallocate(delsq_u)
     deallocate(delsq_circulation_ad)
     deallocate(delsq_circulation)
     deallocate(delsq_vorticity_ad)
     deallocate(delsq_vorticity)
  end if ! (config_h_mom_eddy_visc4 > 0.0)

  if (config_h_mom_eddy_visc2 > 0.0) then
     do iEdge = 1, nEdges
        cell1 = cellsOnEdge(1,iEdge)
        cell2 = cellsOnEdge(2,iEdge)
        vertex1 = verticesOnEdge(1,iEdge)
        vertex2 = verticesOnEdge(2,iEdge)
        do k = 1, nVertLevels
           u_diffusion_ad=tend_u_ad(k,iEdge)
           u_diffusion_ad=meshScalingDel2(iEdge)*config_h_mom_eddy_visc2*u_diffusion_ad
           vorticity_ad(k,vertex1)=vorticity_ad(k,vertex1) +u_diffusion_ad/dvEdge(iEdge)
           vorticity_ad(k,vertex2)=vorticity_ad(k,vertex2) -u_diffusion_ad/dvEdge(iEdge)
           divergence_ad(k,cell1)=divergence_ad(k,cell1) -u_diffusion_ad/dcEdge(iEdge)
           divergence_ad(k,cell2)=divergence_ad(k,cell2) +u_diffusion_ad/dcEdge(iEdge)
           u_diffusion_ad=0.
        end do
     end do ! iEdge = 1, nEdges
  end if ! (config_h_mom_eddy_visc2 > 0.0) then

  do iEdge = nEdges, 1, -1
     cell1 = cellsOnEdge(1,iEdge)
     cell2 = cellsOnEdge(2,iEdge)
     vertex1 = verticesOnEdge(1,iEdge)
     vertex2 = verticesOnEdge(2,iEdge)

     do k = nVertLevels, 1, -1
        h_ad(k,cell1)=h_ad(k,cell1) +gravity/dcEdge(iEdge) * tend_u_ad(k,iEdge)
        h_ad(k,cell2)=h_ad(k,cell2) -gravity/dcEdge(iEdge) * tend_u_ad(k,iEdge)
        ke_ad(k,cell1)=ke_ad(k,cell1) + tend_u_ad(k,iEdge)/dcEdge(iEdge)
        ke_ad(k,cell2)=ke_ad(k,cell2) - tend_u_ad(k,iEdge)/dcEdge(iEdge)
        q_ad=tend_u_ad(k,iEdge)
        tend_u_ad(k,iEdge)=0.

        do j = nEdgesOnEdge(iEdge), 1, -1
           eoe = edgesOnEdge(j,iEdge)
           ! ----- forward calculation -----
           workpv = 0.5 * (pv_edge(k,iEdge) + pv_edge(k,eoe))
           ! ----- end forward calculation -----
           h_edge_ad(k,eoe)=h_edge_ad(k,eoe) + weightsOnEdge(j,iEdge)*u(k,eoe)*workpv*q_ad
           workpv_ad=weightsOnEdge(j,iEdge)*u(k,eoe)*h_edge(k,eoe) * q_ad
           u_ad(k,eoe)=u_ad(k,eoe) + weightsOnEdge(j,iEdge)*workpv*h_edge(k,eoe) * q_ad
           pv_edge_ad(k,eoe)=pv_edge_ad(k,eoe) + 0.5*workpv_ad
           pv_edge_ad(k,iEdge)=pv_edge_ad(k,iEdge) + 0.5*workpv_ad
           workpv_ad=0.
        end do
        q_ad=0.
     end do
  end do
  tend_u_ad(:,:) = 0.0

  do iCell = nCells, 1, -1
     do k = nVertLevels, 1, -1
        tend_h_ad(k,iCell)=tend_h_ad(k,iCell)/areaCell(iCell)
     end do
  end do
  do iEdge = 1, nEdges
     cell1 = cellsOnEdge(1,iEdge)
     cell2 = cellsOnEdge(2,iEdge)
     do k = 1, nVertLevels
        flux_ad=tend_h_ad(k,cell2)
        flux_ad=flux_ad-tend_h_ad(k,cell1)
        h_edge_ad(k,iEdge)=h_edge_ad(k,iEdge) +u(k,iEdge)*dvEdge(iEdge)*flux_ad
        u_ad(k,iEdge)=u_ad(k,iEdge) +dvEdge(iEdge)*h_edge(k,iEdge)*flux_ad
        flux_ad=0.
     end do
  end do
  tend_h_ad(:,:) = 0.0

end subroutine sw_compute_tend_adj

subroutine mpas_reconstruct_2d_adj(nCells, nEdges, maxEdges, nVertLevels, R3, &
     on_a_sphere, edgesOnCell, nEdgesOnCell, latCell, lonCell, &
     coeffs_reconstruct, u, u_ad, uReconstructX, uReconstructX_ad, &
     uReconstructY, uReconstructY_ad, uReconstructZ, uReconstructZ_ad, &
     uReconstructZonal, uReconstructZonal_ad, &
     uReconstructMeridional, uReconstructMeridional_ad)
  implicit none

  integer, parameter :: RKIND  = selected_real_kind(12)

  integer :: nCells, nEdges, nVertLevels, maxEdges, R3
  real (kind=RKIND), dimension(nVertLevels,nEdges), intent(in) :: u
  real (kind=RKIND), dimension(nVertLevels,nEdges) :: u_ad
  real (kind=RKIND), dimension(nVertLevels,nCells), intent(in) :: uReconstructX
  real (kind=RKIND), dimension(nVertLevels,nCells) :: uReconstructX_ad
  real (kind=RKIND), dimension(nVertLevels,nCells), intent(in) :: uReconstructY
  real (kind=RKIND), dimension(nVertLevels,nCells) :: uReconstructY_ad
  real (kind=RKIND), dimension(nVertLevels,nCells), intent(in) :: uReconstructZ
  real (kind=RKIND), dimension(nVertLevels,nCells) :: uReconstructZ_ad
  real (kind=RKIND), dimension(nVertLevels,nCells), intent(in) :: uReconstructZonal
  real (kind=RKIND), dimension(nVertLevels,nCells) :: uReconstructZonal_ad
  real (kind=RKIND), dimension(nVertLevels,nCells), intent(in) :: uReconstructMeridional
  real (kind=RKIND), dimension(nVertLevels,nCells) :: uReconstructMeridional_ad
  !logical, optional, intent(in) :: includeHalos

  integer, dimension(maxEdges,nCells) :: edgesOnCell
  integer, dimension(nCells) :: nEdgesOnCell
  real(kind=RKIND), dimension(nCells) :: latCell, lonCell
  real (kind=RKIND), dimension(R3,maxEdges,nCells) :: coeffs_reconstruct
  logical :: on_a_sphere

  !f2py intent(in) nCells, nEdges, nVertLevels, maxEdges, R3, on_a_sphere,
  !f2py intent(in) edgesOnCell, nEdgesOnCell, latCell, lonCell,
  !f2py intent(in) coeffs_reconstruct, u, uReconstructX, uReconstructY,
  !f2py intent(in) uReconstructZ, uReconstructZonal, uReconstructMeridional

  !f2py intent(in,out) u_ad, uReconstructX_ad, uReconstructY_ad, uReconstructZ_ad,
  !f2py intent(in,out) uReconstructZonal_ad, uReconstructMeridional_ad

  ! local variable
  logical :: includeHalosLocal
  integer :: iCell,iEdge, i
  real (kind=RKIND) :: clat, slat, clon, slon

  if (on_a_sphere) then
     do iCell = nCells, 1, -1
        clat = cos(latCell(iCell))
        slat = sin(latCell(iCell))
        clon = cos(lonCell(iCell))
        slon = sin(lonCell(iCell))

        uReconstructZ_ad(:,iCell)=uReconstructZ_ad(:,iCell) &
             +clat*uReconstructMeridional_ad(:,iCell)
        uReconstructY_ad(:,iCell)=uReconstructY_ad(:,iCell) &
             -slon*slat*uReconstructMeridional_ad(:,iCell)
        uReconstructX_ad(:,iCell)=uReconstructX_ad(:,iCell) &
             -clon*slat*uReconstructMeridional_ad(:,iCell)
        uReconstructMeridional_ad(:,iCell)=0._RKIND

        uReconstructY_ad(:,iCell)=uReconstructY_ad(:,iCell) &
             +clon*uReconstructZonal_ad(:,iCell)
        uReconstructX_ad(:,iCell)=uReconstructX_ad(:,iCell) &
             -slon*uReconstructZonal_ad(:,iCell)
        uReconstructZonal_ad(:,iCell)=0._RKIND
     end do
  else
     do iCell = nCells, 1, -1
        uReconstructY_ad(:,iCell)=uReconstructY_ad(:,iCell) &
             +uReconstructMeridional_ad(:,iCell)
        uReconstructMeridional_ad(:,iCell)=0._RKIND
        uReconstructX_ad(:,iCell)=uReconstructX_ad(:,iCell) &
             +uReconstructZonal_ad(:,iCell)
        uReconstructZonal_ad(:,iCell)=0._RKIND
     end do
  end if

  do iCell = nCells, 1, -1
     do i=nEdgesOnCell(iCell), 1, -1
        iEdge = edgesOnCell(i,iCell)
        u_ad(:,iEdge)=u_ad(:,iEdge) &
             + coeffs_reconstruct(3,i,iCell)*uReconstructZ_ad(:,iCell)
        u_ad(:,iEdge)=u_ad(:,iEdge) &
             + coeffs_reconstruct(2,i,iCell)*uReconstructY_ad(:,iCell)
        u_ad(:,iEdge)=u_ad(:,iEdge) &
             + coeffs_reconstruct(1,i,iCell)*uReconstructX_ad(:,iCell)
     end do
     uReconstructX_ad(:,iCell) = 0.0_RKIND
     uReconstructY_ad(:,iCell) = 0.0_RKIND
     uReconstructZ_ad(:,iCell) = 0.0_RKIND
  end do

end subroutine mpas_reconstruct_2d_adj


text:
*> \brief  SSYSV computes the solution to system of linear equations A * X = B for SY matrices
*
*  =========== DOCUMENTATION ===========
*
* Online html documentation available at
*            http://www.netlib.org/lapack/explore-html/
*
*> \htmlonly
*> Download SSYSV + dependencies
*>
*> [TGZ]
*>
*> [ZIP]
*>
*> [TXT]
*> \endhtmlonly
*
*  Definition:
*  ===========
*
*       SUBROUTINE SSYSV( UPLO, N, NRHS, A, LDA, IPIV, B, LDB, WORK,
*                         LWORK, INFO )
*
*       .. Scalar Arguments ..
*       CHARACTER          UPLO
*       INTEGER            INFO, LDA, LDB, LWORK, N, NRHS
*       ..
*       .. Array Arguments ..
*       INTEGER            IPIV( * )
*       REAL               A( LDA, * ), B( LDB, * ), WORK( * )
*       ..
*
*
*> \par Purpose:
*  =============
*>
*> \verbatim
*>
*> SSYSV computes the solution to a real system of linear equations
*>    A * X = B,
*> where A is an N-by-N symmetric matrix and X and B are N-by-NRHS
*> matrices.
*>
*> The diagonal pivoting method is used to factor A as
*>    A = U * D * U**T,  if UPLO = 'U', or
*>    A = L * D * L**T,  if UPLO = 'L',
*> where U (or L) is a product of permutation and unit upper (lower)
*> triangular matrices, and D is symmetric and block diagonal with
*> 1-by-1 and 2-by-2 diagonal blocks.  The factored form of A is then
*> used to solve the system of equations A * X = B.
*> \endverbatim
*
*  Arguments:
*  ==========
*
*> \param[in] UPLO
*> \verbatim
*>          UPLO is CHARACTER*1
*>          = 'U':  Upper triangle of A is stored;
*>          = 'L':  Lower triangle of A is stored.
*> \endverbatim
*>
*> \param[in] N
*> \verbatim
*>          N is INTEGER
*>          The number of linear equations, i.e., the order of the
*>          matrix A.  N >= 0.
*> \endverbatim
*>
*> \param[in] NRHS
*> \verbatim
*>          NRHS is INTEGER
*>          The number of right hand sides, i.e., the number of columns
*>          of the matrix B.  NRHS >= 0.
*> \endverbatim
*>
*> \param[in,out] A
*> \verbatim
*>          A is REAL array, dimension (LDA,N)
*>          On entry, the symmetric matrix A.  If UPLO = 'U', the leading
*>          N-by-N upper triangular part of A contains the upper
*>          triangular part of the matrix A, and the strictly lower
*>          triangular part of A is not referenced.  If UPLO = 'L', the
*>          leading N-by-N lower triangular part of A contains the lower
*>          triangular part of the matrix A, and the strictly upper
*>          triangular part of A is not referenced.
*>
*>          On exit, if INFO = 0, the block diagonal matrix D and the
*>          multipliers used to obtain the factor U or L from the
*>          factorization A = U*D*U**T or A = L*D*L**T as computed by
*>          SSYTRF.
*> \endverbatim
*>
*> \param[in] LDA
*> \verbatim
*>          LDA is INTEGER
*>          The leading dimension of the array A.  LDA >= max(1,N).
*> \endverbatim
*>
*> \param[out] IPIV
*> \verbatim
*>          IPIV is INTEGER array, dimension (N)
*>          Details of the interchanges and the block structure of D, as
*>          determined by SSYTRF.  If IPIV(k) > 0, then rows and columns
*>          k and IPIV(k) were interchanged, and D(k,k) is a 1-by-1
*>          diagonal block.  If UPLO = 'U' and IPIV(k) = IPIV(k-1) < 0,
*>          then rows and columns k-1 and -IPIV(k) were interchanged and
*>          D(k-1:k,k-1:k) is a 2-by-2 diagonal block.  If UPLO = 'L' and
*>          IPIV(k) = IPIV(k+1) < 0, then rows and columns k+1 and
*>          -IPIV(k) were interchanged and D(k:k+1,k:k+1) is a 2-by-2
*>          diagonal block.
*> \endverbatim
*>
*> \param[in,out] B
*> \verbatim
*>          B is REAL array, dimension (LDB,NRHS)
*>          On entry, the N-by-NRHS right hand side matrix B.
*>          On exit, if INFO = 0, the N-by-NRHS solution matrix X.
*> \endverbatim
*>
*> \param[in] LDB
*> \verbatim
*>          LDB is INTEGER
*>          The leading dimension of the array B.  LDB >= max(1,N).
*> \endverbatim
*>
*> \param[out] WORK
*> \verbatim
*>          WORK is REAL array, dimension (MAX(1,LWORK))
*>          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
*> \endverbatim
*>
*> \param[in] LWORK
*> \verbatim
*>          LWORK is INTEGER
*>          The length of WORK.  LWORK >= 1, and for best performance
*>          LWORK >= max(1,N*NB), where NB is the optimal blocksize for
*>          SSYTRF.
*>          for LWORK < N, TRS will be done with Level BLAS 2
*>          for LWORK >= N, TRS will be done with Level BLAS 3
*>
*>          If LWORK = -1, then a workspace query is assumed; the routine
*>          only calculates the optimal size of the WORK array, returns
*>          this value as the first entry of the WORK array, and no error
*>          message related to LWORK is issued by XERBLA.
*> \endverbatim
*>
*> \param[out] INFO
*> \verbatim
*>          INFO is INTEGER
*>          = 0: successful exit
*>          < 0: if INFO = -i, the i-th argument had an illegal value
*>          > 0: if INFO = i, D(i,i) is exactly zero.  The factorization
*>               has been completed, but the block diagonal matrix D is
*>               exactly singular, so the solution could not be computed.
*> \endverbatim
*
*  Authors:
*  ========
*
*> \author Univ. of Tennessee
*> \author Univ. of California Berkeley
*> \author Univ. of Colorado Denver
*> \author NAG Ltd.
*
*> \ingroup realSYsolve
*
*  =====================================================================
      SUBROUTINE SSYSV( UPLO, N, NRHS, A, LDA, IPIV, B, LDB, WORK,
     $                  LWORK, INFO )
*
*  -- LAPACK driver routine --
*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
*
*     .. Scalar Arguments ..
      CHARACTER          UPLO
      INTEGER            INFO, LDA, LDB, LWORK, N, NRHS
*     ..
*     .. Array Arguments ..
      INTEGER            IPIV( * )
      REAL               A( LDA, * ), B( LDB, * ), WORK( * )
*     ..
*
*  =====================================================================
*
*     .. Local Scalars ..
      LOGICAL            LQUERY
      INTEGER            LWKOPT
*     ..
*     .. External Functions ..
      LOGICAL            LSAME
      EXTERNAL           LSAME
*     ..
*     .. External Subroutines ..
      EXTERNAL           XERBLA, SSYTRF, SSYTRS, SSYTRS2
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC          MAX
*     ..
*     .. Executable Statements ..
*
*     Test the input parameters.
*
      INFO = 0
      LQUERY = ( LWORK.EQ.-1 )
      IF( .NOT.LSAME( UPLO, 'U' ) .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
         INFO = -1
      ELSE IF( N.LT.0 ) THEN
         INFO = -2
      ELSE IF( NRHS.LT.0 ) THEN
         INFO = -3
      ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
         INFO = -5
      ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
         INFO = -8
      ELSE IF( LWORK.LT.1 .AND. .NOT.LQUERY ) THEN
         INFO = -10
      END IF
*
      IF( INFO.EQ.0 ) THEN
         IF( N.EQ.0 ) THEN
            LWKOPT = 1
         ELSE
            CALL SSYTRF( UPLO, N, A, LDA, IPIV, WORK, -1, INFO )
            LWKOPT = WORK(1)
         END IF
         WORK( 1 ) = LWKOPT
      END IF
*
      IF( INFO.NE.0 ) THEN
         CALL XERBLA( 'SSYSV ', -INFO )
         RETURN
      ELSE IF( LQUERY ) THEN
         RETURN
      END IF
*
*     Compute the factorization A = U*D*U**T or A = L*D*L**T.
*
      CALL SSYTRF( UPLO, N, A, LDA, IPIV, WORK, LWORK, INFO )
      IF( INFO.EQ.0 ) THEN
*
*        Solve the system A*X = B, overwriting B with X.
*
         IF ( LWORK.LT.N ) THEN
*
*        Solve with TRS ( Use Level BLAS 2)
*
            CALL SSYTRS( UPLO, N, NRHS, A, LDA, IPIV, B, LDB, INFO )
*
         ELSE
*
*        Solve with TRS2 ( Use Level BLAS 3)
*
            CALL SSYTRS2( UPLO,N,NRHS,A,LDA,IPIV,B,LDB,WORK,INFO )
*
         END IF
*
      END IF
*
      WORK( 1 ) = LWKOPT
*
      RETURN
*
*     End of SSYSV
*
      END

ext: f

max_stars_repo_name: sergey-v-kuznetsov/lapack

max_issues_repo_name: sergey-v-kuznetsov/lapack

max_forks_repo_name: sergey-v-kuznetsov/lapack


text:
      subroutine nplotter_generic(p,wt,wt2,switch)
      implicit none
      include 'vegas_common.f'
      include 'bbproc.f'
      include 'clustering.f'
      include 'constants.f'
      include 'histo.f'
      include 'jetlabel.f'
      include 'npart.f'
      include 'process.f'
      include 'removebr.f'
      include 'nodecay.f'
      include 'useet.f'
      include 'plabel.f'
      include 'outputflags.f'
      include 'nqcdjets.f'
      include 'nproc.f'
cz Add single-top b fraction, Z. Sullivan 1/25/05
      double precision bwgt
      common/btagging/ bwgt

cz //

      integer n,switch,i5,i6,i7,nu,nplotmax
      character tag*4
      double precision DETABB,DPHIZJ
     & ,DPHIBB
     & ,ETAB1
     & ,ETAB2
     & ,ETANOB
     & ,ETARAP
     & ,ETARAPTWO
     & ,ETARAPTHREE
     & ,ETBIN
     & ,ETDOUBLEBIN
     & ,GETET
     & ,M56CLUST
     & ,MBB
     & ,PT
     & ,PTB1
     & ,PTB2
     & ,PTNOB
     & ,PTQ1
     & ,ETAQ1
     & ,PTOTHER
     & ,ETAOTHER
     & ,PTTWO
     & ,PTTHREE
     & ,R
     & ,R57
     & ,R67
     & ,RBB
     & ,SWAP
     & ,WT
     & ,WT2
     & ,YRAP
     & ,YRAPTWO
     & ,YRAPTHREE
     & ,DPHI_LL
     & ,M_LL
     & ,MTRANS
     & ,SCUT1
     & ,SCUT2
     & ,PT34_VETO
     & ,PTQ1_VETO
     & ,ETAQ1_VETO
     & ,PHI56
     & ,PTLEADINGB
     & ,PTLEADINGNONB
c     & ,HT_STOP
c     & ,MERECON
c     & ,MLBNU
c     & ,QETA
c     & ,RECONCORR

cz Added by Z. Sullivan 1/25/05
cz jet(3) will be filled with pt-ordered jets in all processes
cz ibbar  will hold an index to the extra b in t-channel single-top for
cz          use with bwgt
      integer jet(mxpart),jetstart,ibbar,iz,izj,iztmp,
     . inotb,ilight1,ilight2
cz //

      double precision wtbbar,wtnotb,wtlight1,wtlight2
      double precision m34,etmiss,misset,m345678,
     . p(mxpart,4),fphi,HT,etcharm,deltaeta,cosdeltaphi
      double precision eta3,eta4,eta5,eta6,eta7,eta8,eta9,eta10,
     . eta34,eta56,y34
      double precision r45,r56,m345,m348,m567,m678
      double precision pt3,pt4,pt5,pt6,pt7,pt8,pt9,pt10,
     . pt34,pt56,oldpt(5:7)
      double precision pt345,eta345,y345,y3,y4
      double precision bclustmass,etvec(4),tmp5(4),tmp6(4),tmp7(4)
c      double precision langle,plep(4),plep_wrest(4),pw(4),pw_wrest(4)
      integer ssi(4),tmpi,j,k
      double precision ssd(4),tmpd
      integer eventpart,ib1,ib2
      logical first,jetmerge
      logical jetevent
      character*2 ptet
      common/nplotmax/nplotmax
      common/jetmerge/jetmerge
      common/hwwvars/dphi_ll,m_ll,mtrans,scut1,scut2
c      double precision ycut001,ycut005,ycut010,ycut015,ycut020,
c     . y32,y43,z3,z4,z5,z6
      data first/.true./
      save first
ccccc!$omp threadprivate(first,/nplotmax/)
!$omp threadprivate(/jetmerge/)

c--- Check for special plotting instructions
c      if (runstring(1:4) .eq. 'cfmt') then
c        call plots_stop_cfmt(p,wt,wt2)
c      return
c      endif

c--- Set up string for pt or Et
      if (useEt) then
        ptet='Et'
      else
        ptet='pt'
      endif

      if (first) then
        tag='book'
c--- ensure we initialize all possible histograms
        eventpart=npart+3
        eta3=1d3
        pt3=0d0
        eta4=1d3
        pt4=0d0
        eta5=1d3
        pt5=0d0
        eta6=1d3
        pt6=0d0
        eta7=1d3
        pt7=0d0
        eta8=1d3
        pt8=0d0
        eta9=1d3
        pt9=0d0
        eta10=1d3
        pt10=0d0
        eta34=1d3
        y34=1d3
        pt34=0d0
        r45=0d0
        eta56=1d3
        pt56=0d0
        m56clust=0d0
        r56=0d0
        r57=0d0
        r67=0d0
        misset=0d0
        etbin=0d0
        mbb=0d0
        etab1=1d3
        etab2=1d3
        etanob=1d3
        ptb1=0d0
        ptb2=0d0
        ptQ1=0d0
        etaQ1=0d0
        ptother=0d0
        etaother=1d3
        ptnob=0d0
        rbb=0d0
        detabb=1d3
        dphibb=1d3
        m345=0d0
        m348=0d0
        m567=0d0
        m678=0d0
        m345678=0d0
        HT=0d0
      deltaeta=99d0
      cosdeltaphi=99d0
      ssi(1)=5
      ssi(2)=6
        jetmerge=.true.
        jets=nqcdjets
        goto 99
      else
        tag='plot'
      endif

c--- 'eventpart' will contain the number of actual particles that have
c--- a defined momentum. For most processes, this is calculated as follows:
c----  for lowest order and virtual terms switch=0 and eventpart=npart+2
c---   for real events switch=0 and eventpart=npart+2
c---   for real counter-events switch=1 and eventpart=npart+1
c--- There are some processes for which this is not correct and these
c---  are handled with reference to nproc
      eventpart=npart-switch+2

      if (jets .gt. 0) then
        eventpart=4+jets
      endif
      if (    (case .eq. 'W_only') .or. (case .eq. 'Z_only')
     .    .or.(case .eq. 'ggfus0') .or. (case .eq. 'qg_tbq')) then
        eventpart=4+jets
      elseif (case .eq. 'bq_tpq') then
        eventpart=5+jets
      elseif ((case .eq. 'WWqqbr') .or. (case .eq. 'WWnpol')
     .    .or.(case .eq. 'WZbbar') .or. (case .eq. 'ZZlept')
     .    .or.(case .eq. 'HWW_4l') .or. (case .eq. 'HZZ_4l')
     .    .or.(case .eq. 'HWW2lq') .or. (case .eq. 'ggWWbx')
     .    .or.(case .eq. 'HWW_tb') .or. (case .eq. 'HWWint')
     .    .or.(case .eq. 'HWWH+i') .or. (case .eq. 'ggWW4l')
     .    .or.(case .eq. 'HZZ_tb') .or. (case .eq. 'HZZint')
     .    .or.(case .eq. 'HWWjet') .or. (case .eq. 'qq_HWW')
     .    .or.(case .eq. 'WW_jet') .or. (case .eq. 'ZZ_jet')
     .    .or.(case .eq. 'HWWjet') .or. (case .eq. 'HZZjet')
     .    .or.(case .eq. 'HWW2jt') .or. (case .eq. 'HZZ2jt')
     .    .or.(case .eq. 'HWW3jt') .or. (case .eq. 'HZZ3jt')) then
        eventpart=6+jets
      elseif (((case .eq. 'WHbbar') .or. (case .eq. 'ZHbbar'))
     .         .and. (removebr)) then
        eventpart=6+jets
      elseif ((case .eq. 'tt_bbl') .or. (case .eq. 'tt_bbh')
     .   .or. (case .eq. 'tt_bbu') .or. (case .eq. 'tt_ldk')
     .   .or. (case .eq. 'tt_udk')) then
        eventpart=min(8+jets,9)
      elseif (case .eq. 'qq_ttg') then
        eventpart=9
      elseif ((case .eq. 'W_twdk') .or. (case .eq. 'Wtdkay')
     .    .or.(case .eq. 'Wtbwdk')) then
        eventpart=6+jets
        if (removebr) eventpart=eventpart+1
      elseif ((case .eq. 'WH__WW') .or. (case .eq. 'ZH__WW')) then
        eventpart=8+jets
      endif
      if (nproc .eq. 73) eventpart=4+jets
      if (case .eq. 'Hi_Zga') eventpart=5+jets
      if (case .eq. 'qq_ttw') eventpart=11

c--- this variable should be set to .true. when the jets are reordered
c---  according to their pt (or Et)
      jetevent=.false.
c--- re-order jets according to pt, for a W/Z/H+jet event
      if ((case .eq. 'W_1jet') .or. (case .eq. 'Z_1jet') .or.
     .    (case .eq. 'W_2jet') .or. (case .eq. 'Z_2jet') .or.
     .    (case .eq. 'W_3jet') .or. (case .eq. 'Z_3jet') .or.
     .    (case .eq. 'ggfus1') .or. (case .eq. 'ggfus2') .or.
     .    (case .eq. 'ggfus3') .or. (case .eq. 'qq_Hqq') .or.
     .    (case .eq. 'qqHqqg') .or. (case .eq. 'qg_tbq') .or.
     .    (case .eq. 'gagajj')) then
        jetevent=.true.
        if (algorithm .eq. 'cone') then
          if (jets .gt. 0) pt5=getet(p(5,4),p(5,1),p(5,2),p(5,3))
          if (jets .gt. 1) pt6=getet(p(6,4),p(6,1),p(6,2),p(6,3))
          if (jets .gt. 2) pt7=getet(p(7,4),p(7,1),p(7,2),p(7,3))
        else
          if (jets .gt. 0) pt5=pt(5,p)
          if (jets .gt. 1) pt6=pt(6,p)
          if (jets .gt. 2) pt7=pt(7,p)
        endif
        i5=5
        i6=6
        i7=7
      if (jets .gt. 0) oldpt(5)=pt5
      if (jets .gt. 1) oldpt(6)=pt6
      if (jets .gt. 2) oldpt(7)=pt7
c--- sort for 2 jets
        if (jets .eq. 2) then
          if (pt6 .gt. pt5) then
            i5=6
            i6=5
          endif
        endif
c--- sort for 3 jets
        if (jets .eq. 3) then
          if ((pt5 .gt. pt6) .and. (pt5 .gt. pt7)) then
             i5=5
            if (pt6 .gt. pt7) then
              i6=6
              i7=7
            else
              i6=7
              i7=6
            endif
          endif
          if ((pt6 .gt. pt5) .and. (pt6 .gt. pt7)) then
             i5=6
            if (pt5 .gt. pt7) then
              i6=5
              i7=7
            else
              i6=7
              i7=5
            endif
          endif
          if ((pt7 .gt. pt5) .and. (pt7 .gt. pt6)) then
             i5=7
            if (pt5 .gt. pt6) then
              i6=5
              i7=6
            else
              i6=6
              i7=5
            endif
          endif
        endif
c--- perform exchange
        do nu=1,4
          tmp5(nu)=p(i5,nu)
          tmp6(nu)=p(i6,nu)
          tmp7(nu)=p(i7,nu)
        enddo
        do nu=1,4
          p(5,nu)=tmp5(nu)
          p(6,nu)=tmp6(nu)
          p(7,nu)=tmp7(nu)
        enddo
      if (jets .gt. 0) pt5=oldpt(i5)
      if (jets .gt. 1) pt6=oldpt(i6)
      if (jets .gt. 2) pt7=oldpt(i7)
      endif

      m34=dsqrt((p(3,4)+p(4,4))**2-(p(3,1)+p(4,1))**2
     .         -(p(3,2)+p(4,2))**2-(p(3,3)+p(4,3))**2)

      if (bbproc .and. clustering) then
c--- returns zero cluster mass if two b's are in one jet
        m56clust=dsqrt(bclustmass(p))
      else
        m56clust=dsqrt(max(0d0,
     .   (p(5,4)+p(6,4))**2-(p(5,1)+p(6,1))**2
     .  -(p(5,2)+p(6,2))**2-(p(5,3)+p(6,3))**2))
      endif

      eta3=etarap(3,p)
      y3=yrap(3,p)
      pt3=pt(3,p)
      eta4=etarap(4,p)
      y4=yrap(4,p)
      pt4=pt(4,p)
      eta34=etaraptwo(3,4,p)
      y34=yraptwo(3,4,p)
      pt34=pttwo(3,4,p)
      HT=pt3+pt4

      if (eventpart .gt. 4) then
      eta5=etarap(5,p)
      if (jetevent .eqv. .false.) pt5=pt(5,p)
      r45=R(p,4,5)
      m345=dsqrt((p(3,4)+p(4,4)+p(5,4))**2-(p(3,1)+p(4,1)+p(5,1))**2
     .          -(p(3,2)+p(4,2)+p(5,2))**2-(p(3,3)+p(4,3)+p(5,3))**2)
      pt345=ptthree(3,4,5,p)
      eta345=etarapthree(3,4,5,p)
      y345=yrapthree(3,4,5,p)
      HT=HT+pt5
      endif

      if (eventpart .gt. 5) then
      eta6=etarap(6,p)
      if (jetevent .eqv. .false.) pt6=pt(6,p)
      eta56=etaraptwo(5,6,p)
      eta56=etaraptwo(5,6,p)
      pt56=pttwo(5,6,p)
      r56=R(p,5,6)
      phi56=fphi(5,6,p)
      HT=HT+pt6
      endif

      if (eventpart .gt. 6) then
      eta7=etarap(7,p)
      if (jetevent .eqv. .false.) pt7=pt(7,p)
      r57=R(p,5,7)
      r67=R(p,6,7)
      m567=dsqrt((p(5,4)+p(6,4)+p(7,4))**2-(p(5,1)+p(6,1)+p(7,1))**2
     .          -(p(5,2)+p(6,2)+p(7,2))**2-(p(5,3)+p(6,3)+p(7,3))**2)
      HT=HT+pt7
      endif

      if (eventpart .gt. 7) then
      eta8=etarap(8,p)
      pt8=pt(8,p)
      m348=dsqrt((p(3,4)+p(4,4)+p(8,4))**2-(p(3,1)+p(4,1)+p(8,1))**2
     .          -(p(3,2)+p(4,2)+p(8,2))**2-(p(3,3)+p(4,3)+p(8,3))**2)
      m678=dsqrt((p(6,4)+p(7,4)+p(8,4))**2-(p(6,1)+p(7,1)+p(8,1))**2
     .          -(p(6,2)+p(7,2)+p(8,2))**2-(p(6,3)+p(7,3)+p(8,3))**2)
      m345678=dsqrt((p(3,4)+p(4,4)+p(5,4)+p(6,4)+p(7,4)+p(8,4))**2
     .             -(p(3,1)+p(4,1)+p(5,1)+p(6,1)+p(7,1)+p(8,1))**2
     .             -(p(3,2)+p(4,2)+p(5,2)+p(6,2)+p(7,2)+p(8,2))**2
     .             -(p(3,3)+p(4,3)+p(5,3)+p(6,3)+p(7,3)+p(8,3))**2)
      HT=HT+pt8
      endif

      if (eventpart .gt. 8) then
      eta9=etarap(9,p)
      pt9=pt(9,p)
      HT=HT+pt9
      endif

      if (eventpart .gt. 9) then
      eta10=etarap(10,p)
      pt10=pt(10,p)
      HT=HT+pt10
      endif

      misset=etmiss(p,etvec)

c--- set-up variables to catch b's
        if (bbproc) then
          if     (jets .eq. 1) then
            write(6,*) 'Error: bbproc set, but only 1 jet in nplotter.f'
            stop
          elseif (jets .eq. 2) then
            mbb=m56clust
            ptb1=pt5
            ptb2=pt6
            etab1=eta5
            etab2=eta6
            rbb=r56
            if (ptb2 .gt. ptb1) then
              swap=ptb1
              ptb1=ptb2
              ptb2=swap
              swap=etab1
              etab1=etab2
              etab2=swap
            endif
          elseif (jets .eq. 3) then
            call getbs(p,ib1,ib2)
            if     (ib1 .eq. 5) then
              ptb1=pt5
              etab1=eta5
            elseif (ib1 .eq. 6) then
              ptb1=pt6
              etab1=eta6
            elseif (ib1 .eq. 7) then
              ptb1=pt7
              etab1=eta7
            endif
            if     (ib2 .eq. 5) then
              ptb2=pt5
              etab2=eta5
            elseif (ib2 .eq. 6) then
              ptb2=pt6
              etab2=eta6
            elseif (ib2 .eq. 7) then
              ptb2=pt7
              etab2=eta7
            endif
            if (ptb2 .gt. ptb1) then
              swap=ptb1
              ptb1=ptb2
              ptb2=swap
              swap=etab1
              etab1=etab2
              etab2=swap
            endif
            if     (ib1+ib2 .eq. 11) then
              mbb=dsqrt((p(5,4)+p(6,4))**2-(p(5,1)+p(6,1))**2
     .                 -(p(5,2)+p(6,2))**2-(p(5,3)+p(6,3))**2)
              ptnob=pt7
              etanob=eta7
            elseif (ib1+ib2 .eq. 12) then
              mbb=dsqrt((p(5,4)+p(7,4))**2-(p(5,1)+p(7,1))**2
     .                 -(p(5,2)+p(7,2))**2-(p(5,3)+p(7,3))**2)
              ptnob=pt6
              etanob=eta6
            elseif (ib1+ib2 .eq. 13) then
              mbb=dsqrt((p(6,4)+p(7,4))**2-(p(6,1)+p(7,1))**2
     .                 -(p(6,2)+p(7,2))**2-(p(6,3)+p(7,3))**2)
              ptnob=pt5
              etanob=eta5
            endif
            rbb=r(p,ib1,ib2)
          else
            if     (nproc .eq. 146) then
              call getbs(p,ib1,ib2)
              ptb1=pt(ib1,p)
              ptb2=pt(ib2,p)
              etab1=etarap(ib1,p)
              etab2=etarap(ib2,p)
              rbb=r(p,ib1,ib2)
            elseif (nproc .eq. 103) then
c--- set b-quarks to 5 and 6 arbitrarily: if you're interested in
c--- this process, these should be set up properly
              ib1=5
              ib2=6
              ptb1=pt(ib1,p)
              ptb2=pt(ib2,p)
              etab1=etarap(ib1,p)
              etab2=etarap(ib2,p)
              rbb=r(p,ib1,ib2)
            elseif ((case .eq. 'tt_bbh') .or. (case .eq. 'qq_ttw')) then
c--- set b-quarks to 5 and 6: this should have been ensured by jetreorder.f
              ib1=5
              ib2=6
              ptb1=pt(ib1,p)
              ptb2=pt(ib2,p)
              etab1=etarap(ib1,p)
              etab2=etarap(ib2,p)
              rbb=r(p,ib1,ib2)
            else
            write(6,*) 'Unforeseen # of jets and b-quarks in nplotter.f'
            stop
            endif
          endif
          detabb=dabs(etab1-etab2)
          dphibb=dsqrt(dabs(rbb**2-detabb**2))
        endif

c--- special variables for ordering of b-quarks in ttH production
      if (case .eq. 'qq_tth') then
        eventpart=10 ! ensure all histograms are filled
c--- index into pt-ordered list of b-quarks using ssi(n), where n=1,..,4
c--- ordered list is obtained by shuffle-sort
        ssi(1)=5
      ssi(2)=6
      ssi(3)=9
      ssi(4)=10
      do j=1,4
        ssd(j)=pt(ssi(j),p)
      enddo
      do j=4,1,-1
      do k=j-1,1,-1
        if(ssd(j) .gt. ssd(k)) then
          tmpd=ssd(k)
          tmpi=ssi(k)
          ssd(k)=ssd(j)
          ssi(k)=ssi(j)
          ssd(j)=tmpd
          ssi(j)=tmpi
        endif
      enddo
      enddo
c--- example: invariant mass of two highest-pt jets
        mbb=dsqrt(max(0d0,
     .   (p(ssi(1),4)+p(ssi(2),4))**2-(p(ssi(1),1)+p(ssi(2),1))**2
     .  -(p(ssi(1),2)+p(ssi(2),2))**2-(p(ssi(1),3)+p(ssi(2),3))**2))
c      do j=1,4
c      write(6,*) j,ssd(j),ssi(j)
c      enddo
c      pause
      endif

      if     (case .eq. 'gQ__ZQ') then
c--- for the Z+b process, catch the highest pt heavy quark
        call getptQ1(pt5,pt6,eta5,eta6,ptQ1,etaQ1,1)
      elseif (case .eq. 'H_1jet') then
c--- for the H+b process, catch the highest pt heavy quark
        call getptQ1(pt5,pt6,eta5,eta6,ptQ1,etaQ1,1)
        if (jets .gt. 1) then
          ptother=pt5+pt6-ptQ1
          etaother=eta5+eta6-etaQ1
        else
          ptother=-1d0
          etaother=99d0
        endif
        if (jets .eq. 1) then
          pt34_veto=pt34
          ptQ1_veto=ptQ1
          etaQ1_veto=etaQ1
        else
          pt34_veto=-1d0
          ptQ1_veto=-1d0
          etaQ1_veto=99d0
        endif
      else
        ptQ1=-1d0
        etaQ1=99d0
        ptother=-1d0
        etaother=99d0
      endif

c--- for The Z+Q+jet process, make histograms of the leading jet pt
c---  when it is a b-quark and when it is not
      ptleadingb=-1d0
      ptleadingnonb=-1d0
      if (case .eq. 'Z_bjet') then
        if     ((jetlabel(1) .eq. 'bq') .and. (pt5 .gt. pt6)) then
        ptleadingb=pt5
      elseif ((jetlabel(2) .eq. 'bq') .and. (pt6 .gt. pt5)) then
        ptleadingb=pt6
      elseif ((jetlabel(1) .eq. 'pp') .and. (pt5 .gt. pt6)) then
        ptleadingnonb=pt5
      elseif ((jetlabel(2) .eq. 'pp') .and. (pt6 .gt. pt5)) then
        ptleadingnonb=pt6
      endif
      endif

      if (nproc .eq. 61) then
        etbin=etdoublebin(pt4,pt5)
      endif

c--- find largest rapidity difference between the jets
      deltaeta=99d0
      cosdeltaphi=99d0
      if ((jets .eq. 2) .or. (jets .eq. 3)) then
        if     (jets .eq. 2) then
        i5=5
        i6=6
        elseif (jets .eq. 3) then
          if (abs(eta5-eta6).gt.max(abs(eta5-eta7),abs(eta6-eta7))) then
          i5=5
          i6=6
        else
          if (abs(eta5-eta7).gt.abs(eta6-eta7)) then
            i5=5
            i6=7
          else
            i5=6
            i6=7
          endif
        endif
        endif
      deltaeta=abs(etarap(i5,p)-etarap(i6,p))
      cosdeltaphi=(p(i5,1)*p(i6,1)+p(i5,2)*p(i6,2))
     .           /dsqrt((p(i5,1)**2+p(i5,2)**2)*(p(i6,1)**2+p(i6,2)**2))
        if (cosdeltaphi .lt. -0.999999999D0) cosdeltaphi=-1d0
c      cosdeltaphi=abs(cosdeltaphi)
      endif

   99 continue

c--- make plots for H(->WW)+jet paper
c      if (runstring(1:6) .eq. 'hjetww') then
c        call hwwjetplots(eventpart,tag,p,wt,wt2)
c      first=.false.
c      return
c      endif

cz Added by Z. Sullivan, 1/25/05
c--- further amended by JC, 19/8/08
cz jet(i) will hold a pt-ordered index to jet momenta.  jet(i)=0 by default
cz This is designed to work for processes where jet momenta come last,
cz  specifically single-top processes 160-179.  Check process.DAT for momenta
cz  ordering before using in other processes.  (E.g. ttbar, WW, ZZ are more
cz  complicated, and need an index sorting method.)
      do iz = 1,mxpart
         jet(iz) = 0
      enddo
      if (jets.gt.0) then
c Figure out where jets start. Uses same algorithm as genclust_kt.f
c--- pick out jets: note that we search to npart+2-isub, to get the
c--- number of particles right. Note that isub=0 for all calls except
c--- the dipole contributions, where isub=1.
         do iz=3,eventpart
            if ( (plabel(iz) .eq. 'pp') .or. (plabel(iz) .eq. 'pj')
     .           .or.(plabel(iz) .eq. 'bq') .or. (plabel(iz) .eq. 'ba')
     .           .or.(plabel(iz) .eq. 'qj') ) then
               jetstart=iz
               goto 221
            endif
         enddo

 221     do iz=1,jets                        ! loop over jets
            jet(iz)=jetstart + iz-1          ! append jet
c           Compare iz vs. ordered list, and move up if needed
            if (iz.gt.1) then
               do izj = iz,2,-1
                  if (pt(jet(izj),p).gt.pt(jet(izj-1),p)) then
c                    swap izj,izj-1
                     iztmp=jet(izj-1)
                     jet(izj-1)=jet(izj)
                     jet(izj)=iztmp
                  endif
               enddo
            endif
         enddo
      endif
cz // end jet(3) filling

cz Extract index to b~ (b for t~) in t-channel single-top
cz Events should be plotted with weight = wt*bwgt instead of just wt.
      ibbar=0
      inotb=0
      if (((nproc.eq.161).or.(nproc.eq.166)).and.(jets.gt.0)) then
c         if (bwgt.gt.1d-10) then    ! some reasonable minimal threshold
c   There are 2 cases:
c   If merging did not occur, then only if jetlabel(i)=='pp' is it the b~.
            do iz=1,jets
               if     (jetlabel(jet(iz)-jetstart+1).eq.'pp') then
               ibbar=jet(iz)
               elseif (jetlabel(jet(iz)-jetstart+1).eq.'qj') then
               inotb=jet(iz)
               endif
            enddo
c   If merging occurred, then it matters whether the merged jet was already
c    a bq or not.  There is a dirty test for the simple case if removebr=true
            if (jetmerge) then
               if (jetstart.eq.6) ibbar=jet(1)   ! b~ was merged into jet(1)
               if (jetstart.eq.5) then
c          To do this correctly requires changing the clustering routine.
c           Therefore, we are ignoring this case for now.  You should always
c           use removebr = .true. for comparisons to event generators anyway.
               endif
            endif
c         endif  ! end bwgt > 1d-10
      endif
cz // end finding index to b~ in t-channel single-top
cz // Can now plot pt-ordered jets, and b~ fraction (with wt*bwgt)

      if ((case .eq. 'bq_tpq') .and. (removebr)) then
      ilight1=0
      ilight2=0
      wtbbar=0d0
      wtnotb=0d0
      wtlight1=0d0
      wtlight2=0d0
c--- there are three cases
      if     (jets .eq. 1) then
        if     (jetlabel(1) .eq. 'qj') then
c---   1) only one jet, that is definitely not a b
c---           jetlabel = (qj)
        ibbar=0
c         inotb already set above
        wtnotb=wt
      elseif (jetlabel(1) .eq. 'pp') then
c---   2) only one jet, that could be a b
c---           jetlabel = (pp)
c         ibbar already set above
        inotb=ibbar
        wtbbar=wt*bwgt
        wtnotb=wt*(1d0-bwgt)
        endif
c---   3) two jets, one of which could be a b
c---           jetlabel = (qj,pp) OR (pp,qj)
      elseif (jets .eq. 2) then
c       ibbar already set above
c       inotb already set above
        wtbbar=wt*bwgt
      wtnotb=wt*bwgt
        ilight1=jet(1) ! they are ordered according to pt
        ilight2=jet(2) ! they are ordered according to pt
      wtlight1=wt*(1d0-bwgt)
      wtlight2=wtlight1
      elseif (jets .eq. 0) then
c--- nothing to set: all variables=0 and should only occur when notag=1
        continue
      else
        write(6,*) 'Error: there should be 1 or 2 jets, instead ',jets
      stop
      endif
      if (first) then
c--- on first call, use dummies to ensure histograms are initialized
      ibbar=6
      inotb=7
      ilight1=6
      ilight2=7
      endif
      endif

c--- Book and fill ntuple if that option is set, remembering to divide
c--- by # of iterations now that is handled at end for regular histos
      if (creatent .eqv. .true.) then
        call bookfill(tag,p,wt/dfloat(itmx))
c--- REMOVED - to produce normal histograms as well
c        return
      endif

c--- Otherwise, fill the histograms
      n=nextnplot

c--- SPECIAL PLOTS FOR T-CHANNEL SINGLE TOP
      if ((case .eq. 'bq_tpq') .and. (removebr)) then
c--- RAPIDITY FIRST
c--- leading jet that isn't a b
c---    first account for events with only one non b jet (always occurs)
c---    (always occurs, except if notag=1)
      if (inotb .ne. 0) then
      call bookplot(n,tag,'eta notb 1',
     . etarap(inotb,p),wtnotb,wt2,-8d0,8d0,0.5d0,'lin')
      endif
c---    then account for events with two non b jets (only if ilight1>0)
      if (ilight1 .ne. 0) then
      call bookplot(n,tag,'eta notb 1',
     . etarap(ilight1,p),wtlight1,wt2,-8d0,8d0,0.5d0,'lin')
      endif
      n=n+1
c--- subleading jet that isn't a b (only if ilight1>0)
      if (ilight1 .ne. 0) then
      call bookplot(n,tag,'eta notb 2',
     . etarap(ilight2,p),wtlight2,wt2,-8d0,8d0,0.5d0,'lin')
      endif
      n=n+1
c--- b jet (only if ibbar>0)
      if (ibbar .ne. 0) then
      call bookplot(n,tag,'eta b',
     . etarap(ibbar,p),wtbbar,wt2,-8d0,8d0,0.5d0,'lin')
      endif
      n=n+1

c--- PT SECOND
c--- leading jet that isn't a b
c---    first account for events with only one non b jet (always occurs)
c---    (always occurs, except if notag=1)
      if (inotb .ne. 0) then
      call bookplot(n,tag,'pt notb 1',
     . pt(inotb,p),wtnotb,wt2,0d0,200d0,5d0,'log')
      endif
c---    then account for events with two non b jets (only if ilight1>0)
      if (ilight1 .ne. 0) then
      call bookplot(n,tag,'pt notb 1',
     . pt(ilight1,p),wtlight1,wt2,0d0,200d0,5d0,'log')
      endif
      n=n+1
c--- subleading jet that isn't a b (only if ilight1>0)
      if (ilight1 .ne. 0) then
      call bookplot(n,tag,'pt notb 2',
     . pt(ilight2,p),wtlight2,wt2,0d0,200d0,5d0,'log')
      endif
      n=n+1
c--- b jet (only if ibbar>0)
      if (ibbar .ne. 0) then
      call bookplot(n,tag,'pt b',
     . pt(ibbar,p),wtbbar,wt2,0d0,200d0,5d0,'log')
      endif
      n=n+1
      endif

c      if (runstring(1:5) .eq. 'carlo') then
cc--- Special histograms for comparison with C. Oleari's
cc---  e+e- --> QQbg calculation (Durham jet rates)
c        ycut001=0d0
c      ycut005=0d0
c      ycut010=0d0
c      ycut015=0d0
c      ycut020=0d0
c      if (tag .eq. 'plot') then
c        call durhamalg(p,npart-switch,y32,y43,z3,z4,z5,z6)
c      endif
cc      write(6,'(i3,3e16.6)') switch,y43,y32,wt
c      if ((y32 .gt. 0.01d0) .and. (y43 .lt. 0.01d0)) ycut001=1d0
c      if ((y32 .gt. 0.05d0) .and. (y43 .lt. 0.05d0)) ycut005=1d0
c      if ((y32 .gt. 0.10d0) .and. (y43 .lt. 0.10d0)) ycut010=1d0
c      if ((y32 .gt. 0.15d0) .and. (y43 .lt. 0.15d0)) ycut015=1d0
c      if ((y32 .gt. 0.20d0) .and. (y43 .lt. 0.20d0)) ycut020=1d0
c        call bookplot(n,tag,'ycut=0.01',ycut001,wt,wt2,
c     .                0.5d0,1.5d0,1d0,'lin')
c        n=n+1
c        call bookplot(n,tag,'ycut=0.05',ycut005,wt,wt2,
c     .                0.5d0,1.5d0,1d0,'lin')
c        n=n+1
c        call bookplot(n,tag,'ycut=0.10',ycut010,wt,wt2,
c     .                0.5d0,1.5d0,1d0,'lin')
c        n=n+1
c        call bookplot(n,tag,'ycut=0.15',ycut015,wt,wt2,
c     .                0.5d0,1.5d0,1d0,'lin')
c        n=n+1
c        call bookplot(n,tag,'ycut=0.20',ycut020,wt,wt2,
c     .                0.5d0,1.5d0,1d0,'lin')
c        n=n+1
c        call bookplot(n,tag,'y43      ',y43,wt,wt2,
c     .                0d0,1d0,0.05d0,'lin')
c        n=n+1
c        call bookplot(n,tag,'y43>0.01 ',y43,wt,wt2,
c     .                0.01d0,1.01d0,0.05d0,'lin')
c        n=n+1
c        call bookplot(n,tag,'y43>0.05 ',y43,wt,wt2,
c     .                0.05d0,1d0,0.05d0,'lin')
c        n=n+1
c        call bookplot(n,tag,'y43>0.10 ',y43,wt,wt2,
c     .                0.10d0,1d0,0.05d0,'lin')
c        n=n+1
c        call bookplot(n,tag,'y43>0.15 ',y43,wt,wt2,
c     .                0.15d0,1d0,0.05d0,'lin')
c        n=n+1
c        call bookplot(n,tag,'y43>0.20 ',y43,wt,wt2,
c     .                0.20d0,1d0,0.05d0,'lin')
c        n=n+1
c      endif

c      if (runstring(1:4) .eq. 'stop') then
cc--- Special histograms for single top search
c        if (tag .eq. 'plot') then
cc        call stopcuts(p,eventpart,ht,qeta,mlbnu,merecon,reconcorr)
cc these variables should already be in the common-block
c        else
c          ht_stop=-1d0
c          qeta=99d0
c          mlbnu=-1d0
c          merecon=-1d0
c          reconcorr=-99d0
c        endif
c      call bookplot(n,tag,'HT',HT_stop,wt,wt2,100d0,500d0,10d0,'lin')
c      n=n+1
c      call bookplot(n,tag,'qeta',qeta,wt,wt2,-2.8d0,3.2d0,0.4d0,'lin')
c      n=n+1
c      call bookplot(n,tag,'mlbnu',mlbnu,wt,wt2,0d0,300d0,10d0,'lin')
c      n=n+1
c      call bookplot(n,tag,'reconcorr?',reconcorr,
c     . wt,wt2,-1d0,1d0,1d0,'lin')
c      n=n+1
c      call bookplot(n,tag,'MEn recon',merecon,
c     . wt,wt2,0d0,300d0,10d0,'lin')
c      n=n+1
c      endif

c      if (runstring(1:3) .eq. 'hww') then
cc--- Special histograms for H->WW search
c        if (tag .eq. 'plot') then
cc these variables should already be in the common-block
c        else
c          dphi_ll=-1d0
c          m_ll=-1d0
c          mtrans=-1d0
c          scut1=-1d0
c          scut2=-1d0
c        endif
c      call bookplot(n,tag,'dphi_ll',dphi_ll,wt,wt2,
c     .              0d0,3.142d0,0.1571d0,'lin')
cc      write(6,*) switch,dphi_ll,wt,wt2
c      n=n+1
c      call bookplot(n,tag,'m_ll',m_ll,wt,wt2,0d0,200d0,5d0,'lin')
c      n=n+1
c      call bookplot(n,tag,'mtrans',mtrans,wt,wt2,0d0,250d0,5d0,'lin')
c      n=n+1
c      call bookplot(n,tag,'scut1',scut1,wt,wt2,0.5d0,1.5d0,1d0,'lin')
c      n=n+1
c      call bookplot(n,tag,'scut2',scut2,wt,wt2,0.5d0,1.5d0,1d0,'lin')
c      n=n+1
c      endif

c--- calculate lepton 3 angle in W rest frame
c      do nu=1,4
c      pw(nu)=p(3,nu)+p(4,nu)
c      plep(nu)=p(3,nu)
c      pw_wrest(nu)=0d0
c      enddo
c      pw_wrest(4)=dsqrt(pw(4)**2-pw(1)**2-pw(2)**2-pw(3)**2)
c      call boostx(plep,pw,pw_wrest,plep_wrest)
c      langle=plep_wrest(3)/
c     . dsqrt(plep_wrest(1)**2+plep_wrest(2)**2+plep_wrest(3)**2)
cc      langle=dacos(langle)
c      call bookplot(n,tag,'lepton 3',langle,wt,wt2,
c     . -1d0,1.01d0,0.1d0,'lin')
c      n=n+1

c--- calculate lepton 4 angle in W rest frame
c      do nu=1,4
c      pw(nu)=p(3,nu)+p(4,nu)
c      plep(nu)=p(4,nu)
c      pw_wrest(nu)=0d0
c      enddo
c      pw_wrest(4)=dsqrt(pw(4)**2-pw(1)**2-pw(2)**2-pw(3)**2)
c      call boostx(plep,pw,pw_wrest,plep_wrest)
c      langle=plep_wrest(3)/
c     . dsqrt(plep_wrest(1)**2+plep_wrest(2)**2+plep_wrest(3)**2)
cc      langle=dacos(langle)
c      call bookplot(n,tag,'lepton 4',langle,wt,wt2,
c     . -1d0,1.01d0,0.1d0,'lin')
c      n=n+1

c--- special plots for W+c
      if (case .eq. 'W_cjet') then
        etcharm=-1d0
        if     (jets .eq. 1) then
        if ((jetlabel(1) .eq. 'bq') .or. (jetlabel(1) .eq. 'ba')) then
          etcharm=dsqrt(p(5,1)**2+p(5,2)**2)
     .                  *p(5,4)/dsqrt(p(5,1)**2+p(5,2)**2+p(5,3)**2)
          endif
      elseif (jets .eq. 2) then
        if ((jetlabel(1) .eq. 'bq') .or. (jetlabel(1) .eq. 'ba')) then
          etcharm=dsqrt(p(5,1)**2+p(5,2)**2)
     .                  *p(5,4)/dsqrt(p(5,1)**2+p(5,2)**2+p(5,3)**2)
          endif
        if ((jetlabel(2) .eq. 'bq') .or. (jetlabel(2) .eq. 'ba')) then
          etcharm=dsqrt(p(6,1)**2+p(6,2)**2)
     .                  *p(6,4)/dsqrt(p(6,1)**2+p(6,2)**2+p(6,3)**2)
          endif
      endif
      if ((etcharm .lt. 0d0) .and. (tag .eq. 'plot')) then
        write(6,*) 'Error in nplotter: no charm jet for Et plotting'
        stop
      endif
        call bookplot(n,tag,'Et_c',etcharm,wt,wt2,0d0,100d0,10d0,'log')
        n=n+1
      endif

c--- special plots for W+cc~ (sign of quark works for W+ only)
      if (case .eq. 'Wbbmas') then
        etcharm=-1d0
        if     (jets .eq. 1) then
        if (jetlabel(1) .eq. 'ba') then
          etcharm=dsqrt(p(5,1)**2+p(5,2)**2)
     .                  *p(5,4)/dsqrt(p(5,1)**2+p(5,2)**2+p(5,3)**2)
          endif
      elseif (jets .eq. 2) then
        if (jetlabel(1) .eq. 'ba') then
          etcharm=dsqrt(p(5,1)**2+p(5,2)**2)
     .                  *p(5,4)/dsqrt(p(5,1)**2+p(5,2)**2+p(5,3)**2)
          endif
        if (jetlabel(2) .eq. 'ba') then
          etcharm=dsqrt(p(6,1)**2+p(6,2)**2)
     .                  *p(6,4)/dsqrt(p(6,1)**2+p(6,2)**2+p(6,3)**2)
          endif
      endif
c--- Not an error here: only plot the quark correlated with the W
c      if ((etcharm .lt. 0d0) .and. (tag .eq. 'plot')) then
c        write(6,*) 'Error in nplotter: no charm jet for Et plotting'
c        stop
c      endif
        call bookplot(n,tag,'Et_c',etcharm,wt,wt2,0d0,100d0,10d0,'log')
        n=n+1
      endif

c--- special plots for computing ttbar asymmetry

c--- inclusive case
      call bookplot(n,tag,'y(t)-y(tbar)',y3-y4,wt,wt2,
     & -4d0,4d0,4d0,'lin')
      n=n+1

c--- for bins of |y(t)-y(tbar)| of (0,0.5), (0.5,1), (1,1.5)
      call bookplot(n,tag,'y(t)-y(tbar)',y3-y4,wt,wt2,
     & -1.5d0,1.5d0,0.5d0,'lin')
      n=n+1

c--- for bins of |y(t)-y(tbar)| of (1.5,4.5)
      call bookplot(n,tag,'y(t)-y(tbar)',y3-y4,wt,wt2,
     & -4.5d0,4.5d0,3d0,'lin')
      n=n+1

c--- added extra plot here, for the angle analysis of G. Hesketh et al.
      if ((case .eq. 'Z_1jet') .or. (case .eq. 'Z_2jet')) then
      dphizj=atan2(p(3,1)+p(4,1),p(3,2)+p(4,2))-atan2(p(5,1),p(5,2))
      if (dphizj .gt. pi) dphizj=twopi-dphizj
      if (dphizj .lt. -pi) dphizj=twopi+dphizj
      call bookplot(n,tag,'dphi_zj',dphizj,wt,wt2,-3.1416d0,3.1416d0,
     .              1.5708d-1,'log')
      n=n+1
      endif

      call bookplot(n,tag,'HT',HT,wt,wt2,0d0,500d0,20d0,'lin')
      n=n+1
c --- Histograms to monitor exclusive/inclusive cross-sections:
      if ((jets .eq. nqcdjets)  .or. (tag .eq. 'book')) then
      call bookplot(n,tag,'= #LO j ',0.5d0,wt,wt2,0d0,1d0,1d0,'lin')
      endif
      n=n+1
      if ((jets .ge. nqcdjets)  .or. (tag .eq. 'book')) then
      call bookplot(n,tag,'>= #LO j',0.5d0,wt,wt2,0d0,1d0,1d0,'lin')
      endif
      n=n+1
      if (nodecay .eqv. .false.) then
      call bookplot(n,tag,'eta3',eta3,wt,wt2,-4d0,4d0,0.2d0,'lin')
      n=n+1
      call bookplot(n,tag,'y3',y3,wt,wt2,-4d0,4d0,0.2d0,'lin')
      n=n+1
      call bookplot(n,tag,ptet//'3',pt3,wt,wt2,0d0,150d0,5d0,'log')
      n=n+1
      call bookplot(n,tag,'eta4',eta4,wt,wt2,-4d0,4d0,0.2d0,'lin')
      n=n+1
      call bookplot(n,tag,'y4',y4,wt,wt2,-4d0,4d0,0.2d0,'lin')
      n=n+1
      call bookplot(n,tag,ptet//'4',pt4,wt,wt2,25d0,50d0,0.5d0,'lin')
      n=n+1
      call bookplot(n,tag,ptet//'4',pt4,wt,wt2,0d0,150d0,5d0,'log')
      n=n+1
      call bookplot(n,tag,ptet//'4',pt4,wt,wt2,0d0,600d0,20d0,'log')
      n=n+1
      call bookplot(n,tag,'eta34',eta34,wt,wt2,-4d0,4d0,0.4d0,'lin')
      n=n+1
      call bookplot(n,tag,'y34',y34,wt,wt2,-4d0,4d0,0.4d0,'lin')
      n=n+1
      call bookplot(n,tag,ptet//'34',pt34,wt,wt2,0d0,100d0,2d0,'log')
      call ebookplot(n,tag,pt34,wt)
      n=n+1
      call bookplot(n,tag,ptet//'34',pt34,wt,wt2,0d0,40d0,2d0,'log')
      n=n+1
      call bookplot(n,tag,'m34',m34,wt,wt2,0d0,200d0,5d0,'lin')
      n=n+1
      call bookplot(n,tag,'m34',m34,wt,wt2,0d0,500d0,10d0,'log')
      n=n+1
      endif
      call bookplot(n,tag,'misset',misset,wt,wt2,0d0,100d0,2d0,'lin')
      n=n+1

      if (eventpart .gt. 4) then
      call bookplot(n,tag,'eta5',eta5,wt,wt2,-4d0,4d0,0.2d0,'lin')
      n=n+1
      call bookplot(n,tag,'eta5',eta5,wt,wt2,-10d0,10d0,0.5d0,'lin')
      n=n+1
      call bookplot(n,tag,ptet//'5',pt5,wt,wt2,0d0,500d0,10d0,'log')
      n=n+1
c      call bookplot(n,tag,ptet//'5',pt5,wt,wt2,0d0,100d0,2d0,'log')
c      n=n+1
      call bookplot(n,tag,'pt345',pt345,wt,wt2,0d0,200d0,5d0,'log')
      n=n+1
      call bookplot(n,tag,'eta345',eta345,wt,wt2,-4d0,4d0,0.2d0,'lin')
      n=n+1
      call bookplot(n,tag,'y345',y345,wt,wt2,-4d0,4d0,0.2d0,'lin')
      n=n+1
      call bookplot(n,tag,'m345',m345,wt,wt2,0d0,200d0,5d0,'lin')
      n=n+1
      call bookplot(n,tag,'R45',R45,wt,wt2,0d0,5d0,0.1d0,'lin')
      n=n+1

c--- The selected heavy quark jet for Z+b and H+b
c      call bookplot(n,tag,'pt34_veto',
c     . pt34_veto,wt,wt2,0d0,200d0,10d0,'log')
c      n=n+1
c      call bookplot(n,tag,'ptQ_central',ptQ1,wt,wt2,0d0,200d0,5d0,'log')
c      n=n+1
c      call bookplot(n,tag,'ptQ_veto',
c     . ptQ1_veto,wt,wt2,0d0,200d0,5d0,'log')
c      n=n+1
c      call bookplot(n,tag,'etaQ_central',etaQ1,
c     . wt,wt2,-5d0,5d0,0.4d0,'lin')
c      n=n+1
c      call bookplot(n,tag,'etaQ_veto',
c     . etaQ1_veto,wt,wt2,-4d0,4d0,0.1d0,'lin')
c      n=n+1
c--- The other jet for H+b (which is a heavy quark for process 143)
c      call bookplot(n,tag,'pt other',ptother,wt,wt2,0d0,200d0,5d0,'log')
c      n=n+1

      endif

      if (eventpart .gt. 5) then
c      call bookplot(n,tag,'eta other',etaother,
c     . wt,wt2,-5d0,5d0,0.4d0,'lin')
c      n=n+1
      call bookplot(n,tag,'eta6',eta6,wt,wt2,-4d0,4d0,0.2d0,'lin')
      n=n+1
      call bookplot(n,tag,ptet//'6',pt6,wt,wt2,0d0,500d0,10d0,'log')
      n=n+1
c      call bookplot(n,tag,ptet//'6',pt6,wt,wt2,0d0,200d0,2d0,'log')
c      n=n+1
      call bookplot(n,tag,'eta56',eta56,wt,wt2,-4d0,4d0,0.2d0,'lin')
      n=n+1
      call bookplot(n,tag,ptet//'56',pt56,wt,wt2,10d0,250d0,10d0,'log')
      n=n+1
      call bookplot(n,tag,'r56',r56,wt,wt2,0d0,4d0,0.1d0,'lin')
      n=n+1
      call bookplot(n,tag,'phi56',phi56,wt,wt2,0d0,3.14d0,0.314d0,'lin')
      n=n+1
      call bookplot(n,tag,'m56',m56clust,wt,wt2,0d0,200d0,5d0,'log')
      n=n+1
      call bookplot(n,tag,'m56',m56clust,wt,wt2,0d0,2000d0,50d0,'log')
      n=n+1
      call bookplot(n,tag,'eta5-eta6',eta5-eta6,
     . wt,wt2,-6d0,6d0,0.2d0,'lin')
      n=n+1
      endif

c      if (bbproc) then
cc--- Leading b Et, 5 GeV bins from 15 to 200
c      call bookplot(n,tag,'bjet1 Et',ptb1,wt,wt2,15d0,200d0,5d0,'log')
c      n=n+1
cc--- Second b Et, 5 GeV bins from 15 to 200
c      call bookplot(n,tag,'bjet2 Et',ptb2,wt,wt2,15d0,200d0,5d0,'log')
c      n=n+1
cc--- Dijet invariant mass, 2 b-jets, 10 GeV bins from 0 to 250
c      call bookplot(n,tag,'bb mass',mbb,wt,wt2,0d0,350d0,10d0,'log')
c      n=n+1
cc--- Delta_R(b,b), 2 b-jets, bins of 0.2 from 0.35 to 4.95
c      call bookplot(n,tag,'deltaRbb',rbb,
c     . wt,wt2,0.35d0,4.95d0,0.2d0,'lin')
c      n=n+1
c      call bookplot(n,tag,'etab1',etab1,wt,wt2,-4d0,4d0,0.2d0,'lin')
c      n=n+1
c      call bookplot(n,tag,ptet//'b1',ptb1,wt,wt2,0d0,200d0,5d0,'log')
c      n=n+1
c      call bookplot(n,tag,'etab2',etab2,wt,wt2,-4d0,4d0,0.2d0,'lin')
c      n=n+1
c      call bookplot(n,tag,ptet//'b2',ptb2,wt,wt2,0d0,200d0,5d0,'log')
c      n=n+1
c      call bookplot(n,tag,'mbb',mbb,wt,wt2,0d0,200d0,10d0,'log')
c      n=n+1
c      call bookplot(n,tag,'mbb',mbb,wt,wt2,0d0,200d0,5d0,'log')
c      n=n+1
c      call bookplot(n,tag,'rbb',rbb,wt,wt2,0d0,4d0,0.1d0,'lin')
c      n=n+1
c      call bookplot(n,tag,'detabb',detabb,wt,wt2,0d0,4d0,0.1d0,'lin')
c      n=n+1
c      call bookplot(n,tag,'dphibb',dphibb,wt,wt2,0d0,4d0,0.1d0,'lin')
c      n=n+1
c      endif

      if (nqcdjets .ge. 2) then
      call bookplot(n,tag,'delta(eta) >= 2 jets',deltaeta,
     . wt,wt2,0d0,10d0,0.5d0,'lin')
      n=n+1
      endif

      if (eventpart .gt. 6) then

c      if (bbproc) then
cc--- Non b-jet Et, 5 GeV bins from 15 to 200
c      call bookplot(n,tag,'non-b Et',ptnob,wt,wt2,15d0,200d0,5d0,'log')
c      n=n+1
c      call bookplot(n,tag,'etanob',etanob,wt,wt2,-4d0,4d0,0.2d0,'lin')
c      n=n+1
c      call bookplot(n,tag,ptet//'nob',ptnob,wt,wt2,0d0,200d0,5d0,'log')
c      n=n+1
c      endif

      call bookplot(n,tag,'eta7',eta7,wt,wt2,-4d0,4d0,0.2d0,'lin')
      n=n+1
      call bookplot(n,tag,ptet//'7',pt7,wt,wt2,0d0,100d0,5d0,'lin')
      n=n+1
      call bookplot(n,tag,ptet//'7',pt7,wt,wt2,0d0,600d0,20d0,'log')
      n=n+1
      call bookplot(n,tag,'r57',r57,wt,wt2,0d0,4d0,0.1d0,'lin')
      n=n+1
      call bookplot(n,tag,'r67',r67,wt,wt2,0d0,4d0,0.1d0,'lin')
      n=n+1
      call bookplot(n,tag,'m567',m567,wt,wt2,108d0,308d0,2d0,'log')
      n=n+1
      call bookplot(n,tag,'eta diff',eta7-(eta5+eta6)/2d0,
     . wt,wt2,-6d0,6d0,0.2d0,'lin')
      n=n+1
      endif

      if (eventpart .gt. 7) then
      call bookplot(n,tag,'eta8',eta8,wt,wt2,-4d0,4d0,0.2d0,'lin')
      n=n+1
      call bookplot(n,tag,ptet//'8',pt8,wt,wt2,0d0,100d0,5d0,'lin')
      n=n+1
      call bookplot(n,tag,'m348',m348,wt,wt2,108d0,308d0,2d0,'log')
      n=n+1
      call bookplot(n,tag,'m678',m678,wt,wt2,108d0,308d0,2d0,'log')
      n=n+1
      call bookplot(n,tag,'m345678',m345678,wt,wt2,0d0,600d0,10d0,'lin')
      n=n+1
c--- extra plots for azimuthal angle difference between two charged
c--- leptons, in the top decay case
      dphi_ll=(p(4,1)*p(7,1)+p(4,2)*p(7,2))
     .     /dsqrt((p(4,1)**2+p(4,2)**2)*(p(7,1)**2+p(7,2)**2))
      if (dphi_ll .gt. +0.9999999D0) dphi_ll=+1D0
      if (dphi_ll .lt. -0.9999999D0) dphi_ll=-1D0
      dphi_ll=dacos(dphi_ll)
      call bookplot(n,tag,'dphi_ll',dphi_ll,wt,wt2,
     &              0d0,3.142d0,0.1571d0,'lin')
      n=n+1
      if ((m345678 .lt. 400d0) .or. (tag .eq. 'book')) then
      call bookplot(n,tag,'dphi_ll + cut mtt',dphi_ll,wt,wt2,
     &              0d0,3.142d0,0.1571d0,'lin')
      endif
      n=n+1
      endif

      if (eventpart .gt. 8) then
      call bookplot(n,tag,'eta9',eta9,wt,wt2,-4d0,4d0,0.2d0,'lin')
      n=n+1
      call bookplot(n,tag,ptet//'9',pt9,wt,wt2,0d0,100d0,5d0,'lin')
      n=n+1
      endif

      if (eventpart .gt. 9) then
      call bookplot(n,tag,'eta10',eta10,wt,wt2,-4d0,4d0,0.2d0,'lin')
      n=n+1
      call bookplot(n,tag,ptet//'10',pt10,wt,wt2,0d0,100d0,5d0,'lin')
      n=n+1
      endif

c--- special plots for the ttH process
      if (case .eq. 'qq_tth') then
      call bookplot(n,tag,'mb1b2',mbb,wt,wt2,0d0,200d0,5d0,'log')
      n=n+1
      call bookplot(n,tag,'ptb1',pt(ssi(1),p),
     .              wt,wt2,0d0,100d0,5d0,'lin')
      n=n+1
      call bookplot(n,tag,'ptb2',pt(ssi(2),p),
     .              wt,wt2,0d0,100d0,5d0,'lin')
      n=n+1
      endif

c--- additional plot for WW process
      if (nproc .eq. 61) then
      call bookplot(n,tag,'etbin',etbin,wt,wt2,0.5d0,25.5d0,1d0,'lin')
      n=n+1
      endif

c--- special plots for Z+b jet processes
      if (case .eq. 'Z_bjet') then
      call bookplot(n,tag,ptet//' b leading',ptleadingb,
     .               wt,wt2,0d0,500d0,20d0,'log')
      n=n+1
      call bookplot(n,tag,ptet//' non-b leading',ptleadingnonb,
     .               wt,wt2,0d0,500d0,20d0,'log')
      n=n+1
      endif

c--- plots useful for processes with 2 or more jets only
      if (nqcdjets .ge. 2) then
      call bookplot(n,tag,'delta(eta) >= 2 jets',deltaeta,
     . wt,wt2,0d0,10d0,0.5d0,'lin')
      n=n+1
      if ((jets .eq. 2) .or. (tag .eq. 'book')) then
      call bookplot(n,tag,'delta(eta) - 2 jets',deltaeta,
     . wt,wt2,0d0,10d0,0.5d0,'lin')
      endif
      n=n+1
      if ((jets .eq. 2) .or. (tag .eq. 'book')) then
      call bookplot(n,tag,'for phi: delta(eta) - 2 jets',deltaeta,
     . wt*cosdeltaphi,wt2,0d0,10d0,0.5d0,'lin')
      endif
      n=n+1
      if ((jets .eq. 2) .or. (tag .eq. 'book')) then
      call bookplot(n,tag,'cos delta(phi) - 2 jets',cosdeltaphi,
     . wt,wt2,0d0,1d0,0.05d0,'lin')
      endif
      n=n+1
      if ((jets .eq. 3) .or. (tag .eq. 'book')) then
      call bookplot(n,tag,'delta(eta) - 3 jets',deltaeta,
     . wt,wt2,0d0,10d0,0.5d0,'lin')
      endif
      n=n+1
      if ((jets .eq. 3) .or. (tag .eq. 'book')) then
      call bookplot(n,tag,'for phi: delta(eta) - 3 jets',deltaeta,
     . wt*cosdeltaphi,wt2,0d0,10d0,0.5d0,'lin')
      endif
      n=n+1
      if ((jets .eq. 3) .or. (tag .eq. 'book')) then
      call bookplot(n,tag,'cos delta(phi) - 3 jets',cosdeltaphi,
     . wt,wt2,0d0,1d0,0.05d0,'lin')
      endif
      n=n+1
      endif

      n=n-1

c--- ensure the built-in maximum number of histograms is not exceeded
      call checkmaxhisto(n)

c--- set the maximum number of plots, on the first call
      if (first) then
        first=.false.
        nplotmax=n
      endif

cz
      bwgt=0d0  ! for safety
cz //

      return
      end

      subroutine cross(p,i,j,r)
      implicit none
      include 'constants.f'
      integer i,j
      double precision p(mxpart,4),r(3)

      r(1)=p(i,2)*p(j,3)-p(j,2)*p(i,3)
      r(2)=p(i,3)*p(j,1)-p(j,3)*p(i,1)
      r(3)=p(i,1)*p(j,2)-p(j,1)*p(i,2)

      return
      end

      double precision function fphi(n1,n2,p)
      implicit none
      include 'constants.f'
      integer n1,n2
      double precision p(mxpart,4)

      fphi=p(n1,1)*p(n2,1)+p(n1,2)*p(n2,2)
      fphi=fphi/dsqrt(p(n1,1)**2+p(n1,2)**2)
      fphi=fphi/dsqrt(p(n2,1)**2+p(n2,2)**2)
      if     (fphi .gt. +0.9999999D0) then
        fphi=0d0
      elseif (fphi .lt. -0.9999999D0) then
        fphi=pi
      else
        fphi=dacos(fphi)
      endif

      return
      end

      double precision function coslpairet(n1,n2,nm1,nm2,p)
      implicit none
      include 'constants.f'
      integer n1,n2,nm1,nm2,i
      double precision p(mxpart,4),misset(4),pp(4)

      if (nm2 .eq. 0) then
        do i=1,4
          misset(i)=p(nm1,i)
        enddo
      else
        do i=1,4
          misset(i)=p(nm1,i)+p(nm2,i)
        enddo
      endif

      do i=1,4
        pp(i)=p(n1,i)+p(n2,i)
      enddo

      coslpairet=pp(1)*misset(1)+pp(2)*misset(2)
      coslpairet=coslpairet/dsqrt(pp(1)**2+pp(2)**2)
      coslpairet=coslpairet/dsqrt(misset(1)**2+misset(2)**2)

      return
      end

      double precision function deltar(i,j,p)
      implicit none
      include 'constants.f'
      double precision p(mxpart,4),phi1,phi2,etarap,dphi
      integer i,j

      phi1=atan2(p(i,1),p(i,2))
      phi2=atan2(p(j,1),p(j,2))
      dphi=phi1-phi2
      if (dphi .gt. pi) dphi=twopi-dphi
      if (dphi .lt. -pi) dphi=twopi+dphi
      deltar=(etarap(i,p)-etarap(j,p))**2+dphi**2
      deltar=dsqrt(deltar)

      return
      end

      subroutine checkmaxhisto(n)
c--- ensure the built-in maximum number of histograms is not exceeded
      implicit none
      integer n
      include 'histo.f'

      if (n .gt. maxhisto) then
        write(6,*) 'WARNING - TOO MANY HISTOGRAMS!'
        write(6,*) n,' > ',maxhisto,', which is the built-in maximum.'
      write(6,*) 'To use more histograms, change the value of the'
      write(6,*) 'constant MAXHISTO in src/Inc/nplot.f then do:'
      write(6,*)
      write(6,*) ' make clean; make        to recompile from scratch.'
        write(6,*)
      stop
      endif

      return
      end

max_stars_repo_name: tmartini/JHUGen

max_issues_repo_name: tmartini/JHUGen

max_forks_repo_name: tmartini/JHUGen


text:
! Compile and run:  make run_gauss_jordan_test

program gauss_jordan_test
  use gauss_jordan
  implicit none

  real, dimension(3, 3) :: A
  real, dimension(3) :: y, x, x_exp
  real, dimension(size(A, 1) + 1, size(A, 2)) :: Ay

  A(1, :) = (/  0,  2,  1 /)
  A(2, :) = (/  1, -2, -3 /)
  A(3, :) = (/ -1,  1,  2 /)
  y       = (/ -8,  0,  3 /)
  x_exp   = (/ -4, -5,  2 /)

  print *
  print "(a)", "Solving matrix equation Ax = y using Gauss-Jordan..."
  print *

  print "(a)", "A:"
  print "(3f6.1)", transpose(A)
  print *

  print "(a)", "y:"
  print "(3f6.1)", y
  print *

  Ay(1:size(A, 1), :) = transpose(A)
  Ay(size(Ay, 1), :) = y
  print "(a)", "Augmented [A | y]:"
  print "(4f6.1)", Ay
  print *

  Ay = upper_triangular(Ay)
  print "(a)", "Upper triangular form:"
  print "(4f6.1)", Ay
  print *

  Ay = substitution(Ay)
  print "(a)", "Substitution:"
  print "(4f6.1)", Ay
  print *

  x = solve(A, y)
  print "(a)", "Answer:"
  print "(3f6.1)", x
  print *

  print "(a)", "Expected answer:"
  print "(3f6.1)", x_exp
  print *
end program gauss_jordan_test

max_stars_repo_name: YodaEmbedding/experiments

max_issues_repo_name: YodaEmbedding/experiments

max_forks_repo_name: YodaEmbedding/experiments


text:
      module pmep_C
      USE vast_kind_param, ONLY:  double
      integer, dimension(:), allocatable :: ian
      integer :: natom
      real(double), dimension(:,:), allocatable :: co
      end module pmep_C

max_stars_repo_name: openmopac/MOPAC-archive

max_issues_repo_name: openmopac/MOPAC-archive

max_forks_repo_name: openmopac/MOPAC-archive


text:
/*									tab:8
 *
 * ece220_lex.f - lexer for an ECE190 compiler MP (F08 MP5)
 *
 * "Copyright (c) 2006-2018 by Steven S. Lumetta."
 *
 * Permission to use, copy, modify, and distribute this software and its
 * documentation for any purpose, without fee, and without written
 * agreement is hereby granted, provided that the above copyright notice
 * and the following two paragraphs appear in all copies of this software,
 * that the files COPYING and NO_WARRANTY are included verbatim with
 * any distribution, and that the contents of the file README are included
 * verbatim as part of a file named README with any distribution.
 *
 * IN NO EVENT SHALL THE AUTHOR BE LIABLE TO ANY PARTY FOR DIRECT,
 * INDIRECT, SPECIAL, INCIDENTAL, OR CONSEQUENTIAL DAMAGES ARISING OUT
 * OF THE USE OF THIS SOFTWARE AND ITS DOCUMENTATION, EVEN IF THE AUTHOR
 * HAS BEEN ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 *
 * THE AUTHOR SPECIFICALLY DISCLAIMS ANY WARRANTIES, INCLUDING, BUT NOT
 * LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
 * A PARTICULAR PURPOSE.  THE SOFTWARE PROVIDED HEREUNDER IS ON AN "AS IS"
 * BASIS, AND THE AUTHOR NO OBLIGATION TO PROVIDE MAINTENANCE, SUPPORT,
 * UPDATES, ENHANCEMENTS, OR MODIFICATIONS."
 *
 * Author:	    Steve Lumetta
 * Version:	    2
 * Creation Date:   29 September 2008
 * Filename:	    ece220_lex.f
 * History:
 *	SL	2	Tue Apr 03 23:15:33 2018
 *		Changed 190 to 220 for ECE220.
 *	SSL	1	29 September 2008
 *		Based on lexer for toy architecture.
 */

%option noyywrap nounput

%{

#include <ctype.h>
#include <math.h>
#include <stdarg.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <unistd.h>

#include "ece220_parse.tab.h"
#include "ece220_lex.h"

/* not part of ANSI, but in gcclib */
extern char* strdup (const char* s);
extern int fileno (FILE* stream);
extern int vsnprintf (char* s, size_t len, const char* fmt, ...);

ece220_token_t ece220token;

static int16_t parse_int (const char* txt);

static int32_t line_num;

#define YY_NO_INPUT 1

%}

/* operand types */
DECIMAL  [0-9]+
SDECIMAL [-]?{DECIMAL}

IDENT	 [A-Za-z_][A-Za-z_0-9]*

/* operand and white space specification */
SPACE     [ \t]
ENDLINE   \r?\n\r?

/*
 * exclusive lexing states to read comments and to skip input after error
 */
%x ls_comment ls_error

%%
    /* rules for numeric constants */
{SDECIMAL} {ece220token.ival = parse_int (yytext); return C220_KEY_INT;}

    /* rules for keywords */
else	{return C220_KEY_ELSE;}
for	{return C220_KEY_FOR;}
if	{return C220_KEY_IF;}
int	{return C220_KEY_INT_TYPE;}
return	{return C220_KEY_RETURN;}
DEBUG   {return C220_KEY_DEBUG_MARKER;}

    /* multi-character operators */
"&&"	{return C220_KEY_AND;}
"||"	{return C220_KEY_OR;}
"=="	{return C220_KEY_EQ;}
">="	{return C220_KEY_GE;}
"<="	{return C220_KEY_LE;}
"!="	{return C220_KEY_NE;}
"++"    {return C220_KEY_INCR;}
"--"    {return C220_KEY_DECR;}

    /* built-in "library" functions (simpler for students? or for parser?) */
main	{return C220_KEY_MAIN;}
printf	{return C220_KEY_PRINTF;}
rand	{return C220_KEY_RAND;}
scanf	{return C220_KEY_SCANF;}
srand	{return C220_KEY_SRAND;}

    /*
     * rules for variables (must come after
     * keywords to allow 2-character keywords)
     */
{IDENT} {ece220token.name = strdup (yytext); return C220_KEY_IDENT;}

"/*"           {
                   BEGIN (ls_comment);
	       }
<ls_comment>{ENDLINE} {line_num++;}
<ls_comment>"*/"      {BEGIN (0);}
<ls_comment>.         {}

\"[^\"\n\r]*\"	{ece220token.name = strdup (yytext); return C220_KEY_STRING;}
\"[^\"]		{ece220error ("unterminated string");}

    /* eat white space */
{SPACE}+ {}
{ENDLINE} {line_num++;}
. {return *yytext;}

    /* ignore any text after an error */
<ls_error>.         {}
<ls_error>{ENDLINE} {}
<ls_error><<EOF>>   {return -1;}

%%

static int16_t
parse_int (const char* txt)
{
    int32_t num;

    /* we allow only 16-bit values */
    if (1 != sscanf (txt, "%d", &num)) {
        ece220error ("bad numeric constant");
	return 0;
    }
    if (-32768 > num || 32767 < num) {
        ece220error ("numeric constant too large for 16 bits");
	return 0;
    }
    return num;
}

void
ece220error (const char* msg)
{
    fprintf (stderr, "line %3d: %s\n", line_num, msg);
    BEGIN (ls_error);
}

void
ece220errorf (const char *fmt, ...)
{
    char buf[200];
    va_list args;

    /* Check format parameter. */
    if (fmt == NULL)
        return;

    /* Append formatted string. */
    va_start (args, fmt);
    vsnprintf (buf, sizeof (buf), fmt, args);
    va_end (args);

    /* Send it to stderr. */
    fprintf (stderr, "line %3d: %s\n", line_num, buf);

    BEGIN (ls_error);
}

void
ece220_reset_lexer ()
{
    line_num = 1;
}