KERNEL(caldyn_wflux)
  SEQUENCE_C0
    BODY('llm-1,1,-1')
      ! cumulate mass flux convergence from top to bottom
      convm(CELL) = convm(CELL) + convm(UP(CELL))
    END_BLOCK
    EPILOGUE(1)
      dmass_col(HIDX(CELL)) = convm(CELL)
    END_BLOCK
    BODY('2,llm')
      ! Compute vertical mass flux (l=1,llm+1 set to zero at init)
      wflux(CELL) = mass_bl(CELL) * dmass_col(HIDX(CELL)) - convm(CELL) 
    END_BLOCK
  END_BLOCK
  ! make sure wflux is up to date
  BARRIER
END_BLOCK

KERNEL(caldyn_dmass)
  FORALL_CELLS()
    ON_PRIMAL
      convm(CELL) = mass_dbk(CELL) * dmass_col(HIDX(CELL))
    END_BLOCK
  END_BLOCK
END_BLOCK

KERNEL(caldyn_vert)

  DO iq=1,nqdyn
    FORALL_CELLS('2', 'llm')
      ON_PRIMAL
        dtheta_rhodz(CELL,iq) = dtheta_rhodz(CELL,iq) + 0.5*(theta(CELL,iq)+theta(DOWN(CELL),iq))*wflux(CELL)  
      END_BLOCK
    END_BLOCK
    FORALL_CELLS('1', 'llm-1')
      ON_PRIMAL
        dtheta_rhodz(CELL,iq) = dtheta_rhodz(CELL,iq) - 0.5*(theta(CELL,iq)+theta(UP(CELL),iq))*wflux(UP(CELL))  
      END_BLOCK
    END_BLOCK
  END DO

  IF(caldyn_vert_variant == caldyn_vert_cons) THEN
    ! conservative vertical transport of momentum : (F/m)du/deta = 1/m (d/deta(Fu)-u.dF/deta)
    FORALL_CELLS('2','llm')
      ON_EDGES
        wwuu(EDGE) = .25*(wflux(CELL1)+wflux(CELL2))*(u(EDGE)+u(DOWN(EDGE))) ! Fu
      END_BLOCK
    END_BLOCK
    ! make sure wwuu is up to date
    BARRIER

    FORALL_CELLS()
      ON_EDGES
        dFu_deta = wwuu(UP(EDGE))-wwuu(EDGE)                                          ! d/deta (F*u)
	dF_deta  = .5*(wflux(UP(CELL1))+wflux(UP(CELL2))-(wflux(CELL1)+wflux(CELL2))) ! d/deta(F)
        du(EDGE) = du(EDGE) - (dFu_deta-u(EDGE)*dF_deta) / (.5*(rhodz(CELL1)+rhodz(CELL2)))  ! (F/m)du/deta
      END_BLOCK
    END_BLOCK
  ELSE  
    FORALL_CELLS('2','llm')
      ON_EDGES
        wwuu(EDGE) = .5*(wflux(CELL1)+wflux(CELL2))*(u(EDGE)-u(DOWN(EDGE)))
      END_BLOCK
    END_BLOCK

    ! make sure wwuu is up to date
    BARRIER

    FORALL_CELLS()
      ON_EDGES
        du(EDGE) = du(EDGE) - (wwuu(EDGE)+wwuu(UP(EDGE))) / (rhodz(CELL1)+rhodz(CELL2))
      END_BLOCK
    END_BLOCK
  END IF

END_BLOCK

KERNEL(gradient)
  FORALL_CELLS_EXT()
    ON_EDGES
      grad(EDGE) = SIGN*(b(CELL2)-b(CELL1))
    END_BLOCK
  END_BLOCK
END_BLOCK

KERNEL(div)
  FORALL_CELLS_EXT()
    ON_PRIMAL
      div_ij=0.
      FORALL_EDGES
        div_ij = div_ij + SIGN*LE_DE*u(EDGE)
      END_BLOCK
      divu(CELL) = div_ij / AI
    END_BLOCK
  END_BLOCK
END_BLOCK

KERNEL(theta)
  IF(caldyn_eta==eta_mass) THEN ! Compute mass
    ! FIXME : here mass_col is computed from rhodz
    ! so that the DOFs are the same whatever caldyn_eta
    ! in DYNAMICO mass_col is prognosed rather than rhodz
    SEQUENCE_C1
      PROLOGUE(0)
        mass_col(HIDX(CELL))=0.
      END_BLOCK
      BODY('1,llm')
          mass_col(HIDX(CELL)) = mass_col(HIDX(CELL)) + rhodz(CELL)
      END_BLOCK
    END_BLOCK
    FORALL_CELLS_EXT()
      ON_PRIMAL
        ! FIXME : formula below (used in DYNAMICO) is for dak, dbk based on pressure rather than mass
        !        m = mass_dak(CELL)+(mass_col(HIDX(CELL))*g+ptop)*mass_dbk(CELL)
      	!        rhodz(CELL) = m/g
        rhodz(CELL) = mass_dak(CELL) + mass_col(HIDX(CELL))*mass_dbk(CELL)
      END_BLOCK
    END_BLOCK
  END IF
  DO iq=1,nqdyn
    FORALL_CELLS_EXT()
      ON_PRIMAL
        theta(CELL,iq) = theta_rhodz(CELL,iq)/rhodz(CELL)
      END_BLOCK
    END_BLOCK
  END DO
END_BLOCK