KERNEL(advect_vert)

  FORALL_CELLS_EXT('2','llm')
    ON_PRIMAL
      dzqw(CELL) = q(CELL)-q(DOWN(CELL))
    END_BLOCK
  END_BLOCK

  ! We need a barrier here because we compute dzqw above and do a vertical average below
  BARRIER  

  FORALL_CELLS_EXT()
    ON_PRIMAL
      CST_IFTHEN(IS_BOTTOM_LEVEL || IS_TOP_LAYER)
        dzq(CELL)=0.
      CST_ELSE
        dzq_down=dzqw(CELL)
        dzq_up=dzqw(UP(CELL)
        IF(dzq_up*dzq_down>0.) THEN
           dzqm = .5*(dzqup+dzq_down)
           dzqmax = max_slope * min( ABS(dzqw_down), ABS(dzq_up) )
           dzq(CELL) = sign( min(abs(dzqm),dzqmax) , dzqm )  ! sign(a,b)=a*sign(b) 
        ELSE
          dzq(CELL)=0.
        END IF
      CST_ENDIF
    END_BLOCK
  END_BLOCK
   
  ! We need a barrier here because we compute dzq above and do a vertical average below
  BARRIER

  FORALL_CELLS_EXT('1','llm+1')
    ON_PRIMAL
      CST_IFTHEN(IS_BOTTOM_LEVEL || IS_TOP_INTERFACE)
        wq(CELL)=0.
      CST_ELSE
        w = fac*wfluxt(CELL)
	IF(w>0.) THEN  ! upward transport, upwind side is at level l
	  sigw = w/mass(DOWN(CELL))
	  qq   = q(DOWN(CELL) + 0.5*(1.-sigw)*dzq(DOWN(CELL)) ! qq = q if sigw=1 , qq = q+dzq/2 if sigw=0
        ELSE           ! downward transport, upwind side is at level l+1
          sigw = w/mass(CELL)
          qq   = q(CELL)      -0.5*(1.+sigw)*dzq(ij,l) ! qq = q if sigw=-1 , qq = q-dzq/2 if sigw=0               
	ENDIF
	wq(ij,l) = w*qq																 
      CST_ENDIF
    END_BLOCK
  END_BLOCK

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

  FORALL_CELLS_EXT()
    ON_PRIMAL
      newmass = mass(CELL) + fac*(wflux(CELL)-wflux(UP(CELL)))
      q(CELL) = ( q(CELL)*mass(CELL) + wq(CELL)-wq(UP(CELL)) ) / newmass
      IF(update_mass) mass(CELL) = newmass
    END_BLOCK
  END_BLOCK

END_BLOCK

KERNEL(advect_horiz)
  ! evaluate tracer field at point cc using piecewise linear reconstruction
  ! q(cc)= q0 + gradq.(cc-xyz_i)  with xi centroid of hexagon
  ! sign*hfluxt>0 iff outgoing
  FORALL_CELLS_EXT()
    ON_EDGES
      IF(SIGN*hfluxt(EDGE)>0.) THEN
        qe = qi(CELL1)
        qe = qe + (cc(EDGE,1)-xyz_i(HIDX(CELL1),1))*gradq3d(CELL1,1) 
        qe = qe + (cc(EDGE,2)-xyz_i(HIDX(CELL1),2))*gradq3d(CELL1,2) 
        qe = qe + (cc(EDGE,3)-xyz_i(HIDX(CELL1),3))*gradq3d(CELL1,3) 
      ELSE
        qe = qi(CELL2)
        qe = qe + (cc(EDGE,1)-xyz_i(HIDX(CELL2),1))*gradq3d(CELL2,1) 
        qe = qe + (cc(EDGE,2)-xyz_i(HIDX(CELL2),2))*gradq3d(CELL2,2) 
        qe = qe + (cc(EDGE,3)-xyz_i(HIDX(CELL2),3))*gradq3d(CELL2,3) 
      END IF
      qflux(EDGE) = hfluxt(EDGE)*qe
      IF(diagflux_on) qfluxt(EDGE) = qfluxt(EDGE)+qflux(EDGE)
    END_BLOCK
  END_BLOCK

  ! update q and, if update_mass, update mass
  FORALL_CELLS()
    ON_PRIMAL
      dmass=0.
      dq=0.
      FORALL_EDGES
        dmass = dmass + SIGN*hfluxt(EDGE)
        dq = dq + SIGN*qflux(EDGE)
      END_BLOCK
      newmass = mass(CELL) - dmass/AI
      qi(CELL) = (qi(CELL)*mass(CELL)-dq/AI) / newmass
      IF(update_mass) mass(CELL)=newmass
    END_BLOCK
  END_BLOCK

END_BLOCK

{% macro perot_reconstruction() %}
  ! Perot reconstruction based on Gauss theorem
  ! u = sum( u.edge_normal * edge_length * (edge_midpoint-cell_centroid) ) /cell_area
  FORALL_CELLS()
    ON_PRIMAL
      cx=centroid(HIDX(CELL),1)
      cy=centroid(HIDX(CELL),2)
      cz=centroid(HIDX(CELL),3)
      ux=0. ; uy=0. ; uz=0.
      FORALL_EDGES
        {{ caller() }}
	ux = ux + ue_le*(.5*(xyz_v(HIDX(VERTEX1),1)+xyz_v(HIDX(VERTEX2),1))-cx)
	uy = uy + ue_le*(.5*(xyz_v(HIDX(VERTEX1),2)+xyz_v(HIDX(VERTEX2),2))-cy)
	uz = uz + ue_le*(.5*(xyz_v(HIDX(VERTEX1),3)+xyz_v(HIDX(VERTEX2),3))-cz)
      END_BLOCK
      fac = scale*(1./AI)
      ucenter(CELL,1)=ux*fac
      ucenter(CELL,2)=uy*fac
      ucenter(CELL,3)=uz*fac
    END_BLOCK
  END_BLOCK
{% endmacro %}

KERNEL(wind_centered)
{% call() perot_reconstruction() %}
ue_le = SIGN*ue(EDGE)*le(HIDX(EDGE))
{% endcall %}
END_BLOCK

KERNEL(flux_centered)
! NB : here the input data is a flux and has already the factor l_e in it
! Input data is rescaled by factor scale
{% call() perot_reconstruction() %}
ue_le = SIGN*ue(EDGE)
{% endcall %}
END_BLOCK

/*------------------------------------- Energetics --------------------------------*/

{% macro energy_flux(energy) %}
  FORALL_CELLS_EXT()
    ON_PRIMAL
      {{ caller() }}
      {{ energy }}(CELL) = {{ energy }}(CELL) + frac*rhodz(CELL)*energy
      pk(CELL) = energy
    END_BLOCK
  END_BLOCK
  FORALL_CELLS_EXT()
    ON_EDGES
      {{ energy }}_flux(EDGE) = {{ energy }}_flux(EDGE) + .5*massflux(EDGE)*(pk(CELL1)+pk(CELL2))
    END_BLOCK
  END_BLOCK
{% endmacro %}

KERNEL(energy_fluxes)
  ! First diagnose geopotential and temperature, column-wise
  BARRIER
  SEQUENCE_C1
    PROLOGUE(llm)
      pk(CELL) = ptop + .5*g*rhodz(CELL)
    END_BLOCK
    BODY('llm-1,1,-1')
      pk(CELL) = pk(UP(CELL)) + (.5*g)*( rhodz(CELL)+rhodz(UP(CELL)) )
    END_BLOCK
  END_BLOCK
  ! NB : at this point pressure is stored in array pk
  !      pk then serves as buffer to store temperature
  SELECT CASE(caldyn_thermo)
  CASE(thermo_theta)
    GEOPOT('temp_ik')
      theta_ik = theta_rhodz(CELL,1)/rhodz(CELL)
      theta(CELL) = theta_ik
      temp_ik = theta_ik*(p_ik/preff)**kappa
      gv = (g*Rd)*temp_ik/p_ik
    END_GEOPOT
  CASE(thermo_entropy)
    GEOPOT('temp_ik')
      theta_ik = theta_rhodz(CELL,1)/rhodz(CELL)
      temp_ik = Treff*exp((theta_ik + Rd*log(p_ik/preff))/cpp)
      theta(CELL) = Treff*exp(theta_ik/cpp)
      gv = (g*Rd)*temp_ik/p_ik   ! specific volume v = Rd*T/p
    END_GEOPOT
  END SELECT
  BARRIER
  ! Now accumulate energies and energy fluxes
  ! NB : at this point temperature is stored in array pk
  !      pk then serves as buffer to store energy
  ! enthalpy
  {% call energy_flux('enthalpy') %}
    energy = cpp*pk(CELL)
  {% endcall %}
  ! potential energy
  {% call energy_flux('epot') %}
    energy = .5*(geopot(UP(CELL))+geopot(CELL))
  {% endcall %}
  ! theta
  {% call energy_flux('thetat') %}
    energy = theta(CELL)
  {% endcall %}
  ! kinetic energy
  {% call energy_flux('ekin') %}
    energy=0.d0
    FORALL_EDGES
      energy = energy + le(HIDX(EDGE))*de(HIDX(EDGE))*ue(EDGE)**2
    END_BLOCK
    energy = energy * (.25/AI)
  {% endcall %}
  ! ulon
  {% call energy_flux('ulon') %}
    cx=centroid(HIDX(CELL),1)
    cy=centroid(HIDX(CELL),2)
    cz=centroid(HIDX(CELL),3)
    ux=0. ; uy=0. ; uz=0.
    FORALL_EDGES
      ue_le = SIGN*ue(EDGE)*le(HIDX(EDGE))
      ux = ux + ue_le*(.5*(xyz_v(HIDX(VERTEX1),1)+xyz_v(HIDX(VERTEX2),1))-cx)
      uy = uy + ue_le*(.5*(xyz_v(HIDX(VERTEX1),2)+xyz_v(HIDX(VERTEX2),2))-cy)
      uz = uz + ue_le*(.5*(xyz_v(HIDX(VERTEX1),3)+xyz_v(HIDX(VERTEX2),3))-cz)
    END_BLOCK
    ulon_i = ux*elon_i(HIDX(CELL),1) + uy*elon_i(HIDX(CELL),2) + uz*elon_i(HIDX(CELL),3)
    energy = ulon_i*(1./AI)
  {% endcall %}

END_BLOCK

KERNEL(remap_eta)
! IN  : rhodz, mass_dak, mass_dbk
! TMP : mass_col, cur_lev, new_rhodz_cum
! OUT : rhodz, old_rhodz, eta
    SEQUENCE_C1
      PROLOGUE('1')
        rhodz_cum(CELL)=0.
        cur_lev(HIDX(CELL))=1
        eta(CELL)=1.
        new_rhodz_cum(HIDX(CELL))=0.
      END_BLOCK
      BODY('1,llm')
        rhodz_cum(UP(CELL)) = rhodz_cum(CELL) + rhodz(CELL)
      END_BLOCK
      EPILOGUE('llm+1')
        mass_col(HIDX(CELL)) = rhodz_cum(CELL)
      END_BLOCK

      BODY('1,llm')
        old_rhodz(CELL) = rhodz(CELL)
        rhodz(CELL) = mass_dak(CELL) + mass_col(HIDX(CELL))*mass_dbk(CELL)
        rhodz_cum_target = new_rhodz_cum(HIDX(CELL)) + rhodz(CELL)
        DO level = cur_lev(HIDX(CELL)), llm
           rhodz_cum_levp1 = rhodz_cum(AT_LEVEL(CELL,level+1))
           IF(rhodz_cum_target<=rhodz_cum_levp1) EXIT
        END DO
        IF(level>llm) level=llm
        rhodz_cum_lev = rhodz_cum(AT_LEVEL(CELL,level))
        ! now rhodz_cum_lev <= rhodz_cum_target <= rhodz_cum_levp1  
        cur_lev(HIDX(CELL)) = level
        new_rhodz_cum(HIDX(CELL)) = rhodz_cum_target
        eta(UP(CELL)) = level + (rhodz_cum_target-rhodz_cum_lev)/(rhodz_cum_levp1-rhodz_cum_lev)
      END_BLOCK
    END_BLOCK
END_BLOCK

KERNEL(remap_theta)
! IN  : thetarhodz, eta
! TMP : thetarhodz_cum, new_thetarhodz_cum
! OUT : thetarhodz
    SEQUENCE_C1
      PROLOGUE('1')
        thetarhodz_cum(CELL)=0.
      END_BLOCK
      BODY('1,llm')
        thetarhodz_cum(UP(CELL)) = thetarhodz_cum(CELL) + thetarhodz(CELL)
      END_BLOCK
        
      BODY('1,llm+1')
        X = eta(CELL)
        level = MIN(llm,FLOOR(X)) ! eta=llm+1 => level=llm, X=1
        X = X-level
        new_thetarhodz_cum(CELL) = thetarhodz_cum(AT_LEVEL(CELL,level))+X*thetarhodz(AT_LEVEL(CELL,level))
      END_BLOCK
      BODY('1,llm')
        thetarhodz(CELL) = new_thetarhodz_cum(UP(CELL)) - new_thetarhodz_cum(CELL)
      END_BLOCK
    END_BLOCK
END_BLOCK

KERNEL(remap_u)
! IN  : u, old_rhodz, rhodz, eta
! TMP : urhodz_cum, new_urhodz_cum
! OUT : u
    SEQUENCE_E0
      PROLOGUE('1')
        urhodz_cum(EDGE)=0.
      END_BLOCK
      BODY('1,llm')
        urhodz(EDGE) = u(EDGE)*(old_rhodz(CELL1)+old_rhodz(CELL2))
        urhodz_cum(UP(EDGE)) = urhodz_cum(EDGE) + urhodz(EDGE)
      END_BLOCK
        
      BODY('1,llm+1')
        X = .5*(eta(CELL1)+eta(CELL2))
        level = MIN(llm,FLOOR(X))
        X = X-level
        new_urhodz_cum(EDGE) = urhodz_cum(AT_LEVEL(EDGE,level))+X*urhodz(AT_LEVEL(EDGE,level))
      END_BLOCK
      BODY('1,llm')
        u(EDGE) = (new_urhodz_cum(UP(EDGE)) - new_urhodz_cum(EDGE)) / (rhodz(CELL1)+rhodz(CELL2))
      END_BLOCK
    END_BLOCK
END_BLOCK

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

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

END_BLOCK

KERNEL(curl_laplacian)
  FORALL_CELLS_EXT()
    ON_DUAL
      etav = 0.d0
      FORALL_EDGES
         etav = etav + SIGN*u(EDGE)
      END_BLOCK
      qv(DUAL_CELL) = etav/AV
    END_BLOCK
  END_BLOCK

  FORALL_CELLS_EXT()
    ON_EDGES
       curlcurl(EDGE) = SIGN*(qv(VERTEX1)-qv(VERTEX2))*(1./LE_DE)
    END_BLOCK
  END_BLOCK

END_BLOCK