{% macro compute_p_and_Aik() %}
{% set p_and_c2_from_theta=caller %}
    SEQUENCE_C1
      BODY('1,llm')
          rho_ij = (g*m_ik(CELL))/(Phi_il(UP(CELL))-Phi_il(CELL))
          {{ p_and_c2_from_theta() }}
          A_ik(CELL) = c2_mik*(tau/g*rho_ij)**2
      END_BLOCK
    END_BLOCK
{%- endmacro %}

KERNEL(caldyn_mil)
  FORALL_CELLS_EXT('1', 'llm+1')
    ON_PRIMAL
       CST_IF(IS_BOTTOM_LEVEL,    m_il(CELL) = .5*rhodz(CELL)                     )
       CST_IF(IS_INNER_INTERFACE, m_il(CELL) = .5*(rhodz(CELL)+rhodz(DOWN(CELL))) )
       CST_IF(IS_TOP_INTERFACE,   m_il(CELL) = .5*rhodz(DOWN(CELL))               )
    END_BLOCK
  END_BLOCK
END_BLOCK

KERNEL(caldyn_vert_NH)
  FORALL_CELLS_EXT()
    ON_PRIMAL
      CST_IF(IS_BOTTOM_LEVEL,     w_ij = ( W(CELL)+.5*W(UP(CELL)) )/mass(CELL) )
      CST_IF(IS_INNER_LAYER , w_ij = .5*( W(CELL)+W(UP(CELL)) )/mass(CELL) )
      CST_IF(IS_TOP_LAYER,        w_ij = ( .5*W(CELL)+W(UP(CELL)) )/mass(CELL) )
      wflux_ij = .5*(wflux(CELL)+wflux(UP(CELL)))
      W_etadot(CELL) = wflux_ij*w_ij
      eta_dot(CELL) = wflux_ij / mass(CELL)
      wcov(CELL) = w_ij*(geopot(UP(CELL))-geopot(CELL))
    END_BLOCK
  END_BLOCK
  ! add NH term to du
  FORALL_CELLS()
    ON_EDGES
      du(EDGE) = du(EDGE) - .5*(wcov(CELL1)+wcov(CELL2))*SIGN*(eta_dot(CELL2)-eta_dot(CELL1))
    END_BLOCK
  END_BLOCK
  ! add NH terms to dW, dPhi
  ! FIXME : TODO top and bottom
  FORALL_CELLS('2','llm')
    ON_PRIMAL
      dPhi(CELL)=dPhi(CELL)-wflux(CELL)*(geopot(UP(CELL))-geopot(DOWN(CELL)))/(mass(DOWN(CELL))+mass(CELL))
    END_BLOCK
  END_BLOCK

  ! We need a barrier here because we compute W_etadot above and do a vertical difference below
  BARRIER

  FORALL_CELLS('1','llm+1')
    ON_PRIMAL
      CST_IFTHEN(IS_BOTTOM_LEVEL)
         dW(CELL) = dW(CELL) - W_etadot(CELL)
      CST_ELSEIF(IS_TOP_INTERFACE)
         dW(CELL) = dW(CELL) + W_etadot(DOWN(CELL))
      CST_ELSE
         dW(CELL) = dW(CELL) + W_etadot(DOWN(CELL)) - W_etadot(CELL)
      CST_ENDIF
    END_BLOCK
  END_BLOCK
END_BLOCK