source: codes/icosagcm/devel/Python/src/kernels_transport.jin @ 615

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_EXT
    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