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

KERNEL(caldyn_fast)
!
  SELECT CASE(caldyn_thermo)
  CASE(thermo_boussinesq)
    FORALL_CELLS()
      ON_PRIMAL
        berni(CELL) = pk(CELL)
        ! from now on pk contains the vertically-averaged geopotential
        pk(CELL) = .5*(geopot(CELL)+geopot(UP(CELL)))
      END_BLOCK
    END_BLOCK
  CASE(thermo_theta)
    FORALL_CELLS()
      ON_PRIMAL
        berni(CELL) = .5*(geopot(CELL)+geopot(UP(CELL)))
      END_BLOCK
    END_BLOCK
  CASE(thermo_entropy)
    FORALL_CELLS()
      ON_PRIMAL
        berni(CELL) = .5*(geopot(CELL)+geopot(UP(CELL)))
        berni(CELL) = berni(CELL) + pk(CELL)*(cpp-theta(CELL,1)) ! Gibbs = Cp.T-Ts = T(Cp-s)
      END_BLOCK
    END_BLOCK
  CASE(thermo_variable_Cp)
    ! thermodynamics with variable Cp
    !           Cp(T) = Cp0 * (T/T0)^nu
    ! =>            h = Cp(T).T/(nu+1)

    FORALL_CELLS()
      ON_PRIMAL
        berni(CELL) = .5*(geopot(CELL)+geopot(UP(CELL)))
        cp_ik = cpp*(pk(CELL)/Treff)**nu
        berni(CELL) = berni(CELL) + pk(CELL)*(cp_ik/(nu+1.)-theta(CELL,1)) ! Gibbs = h-Ts = T(Cp/(nu+1)-s)
      END_BLOCK
    END_BLOCK
  CASE DEFAULT
    PRINT *, 'Unsupported value of caldyn_thermo : ',caldyn_thermo  ! FIXME
    STOP
  END SELECT
!
  FORALL_CELLS()
    ON_EDGES
      due = .5*(theta(CELL1,1)+theta(CELL2,1))*(pk(CELL2)-pk(CELL1)) + berni(CELL2)-berni(CELL1)
      du(EDGE) = du(EDGE) - SIGN*due
      u(EDGE) = u(EDGE) + tau*du(EDGE)
    END_BLOCK
  END_BLOCK  
!
END_BLOCK

KERNEL(pvort_only)
  FORALL_CELLS_EXT()
    ON_DUAL
      etav = 0.d0
      FORALL_EDGES
         etav = etav + SIGN*u(EDGE)
      END_BLOCK
      hv=0.
      FORALL_VERTICES
        hv = hv + RIV2*rhodz(VERTEX)
      END_BLOCK
      qv(DUAL_CELL) = (etav + FV*AV )/(hv*AV)
    END_BLOCK
  END_BLOCK

  FORALL_CELLS()
    ON_EDGES
      qu(EDGE)=0.5d0*(qv(VERTEX1)+qv(VERTEX2))
    END_BLOCK
  END_BLOCK
END_BLOCK

KERNEL(coriolis)
!
  DO iq=1,nqdyn
    FORALL_CELLS_EXT()
      ON_EDGES
        Ftheta(EDGE) = .5*(theta(CELL1,iq)+theta(CELL2,iq))*hflux(EDGE)
      END_BLOCK
    END_BLOCK
    FORALL_CELLS()
      ON_PRIMAL
        divF=0.
	FORALL_EDGES
	  divF = divF + Ftheta(EDGE)*SIGN
	END_BLOCK
        dtheta_rhodz(CELL,iq) = -divF / AI
      END_BLOCK
    END_BLOCK
  END DO ! iq
!
  FORALL_CELLS()
    ON_PRIMAL
      divF=0.
      FORALL_EDGES
        divF = divF + hflux(EDGE)*SIGN
      END_BLOCK
      convm(CELL) = -divF / AI
    END_BLOCK
  END_BLOCK
!
  FORALL_CELLS() 
    ON_EDGES
      du_trisk=0.
      FORALL_TRISK
        du_trisk = du_trisk + WEE*hflux(EDGE_TRISK)*(qu(EDGE)+qu(EDGE_TRISK))
      END_BLOCK
      du(EDGE) = du(EDGE) + .5*du_trisk
    END_BLOCK
  END_BLOCK

END_BLOCK