#include "ERF_Diffusion.H"
#include "ERF_EddyViscosity.H"
#include "ERF_SolveTridiag.H"
#include "ERF_GetRhoAlpha.H"
#include "ERF_GetRhoAlphaForFaces.H"
#include "ERF_SetupVertDiff.H"

Include dependency graph for ERF_ImplicitDiff_T.cpp:

Macros
#define	INSTANTIATE_IMPLICIT_DIFF_FOR_MOM_LU(STAGDIR)

Functions
void	ImplicitDiffForStateLU_T (const Box &bx, const Box &domain, const int level, const int n, const double dt_d, const GpuArray< Real, AMREX_SPACEDIM 2 > &bc_neumann_vals, const Array4< Real > &cell_data, const Array4< const Real > &z_nd, const Array4< const Real > &detJ, const GpuArray< Real, AMREX_SPACEDIM > &cellSizeInv, const Array4< const Real > &scalar_zflux, const Array4< const Real > &mu_turb, const SolverChoice &solverChoice, const BCRec bc_ptr, const bool use_SurfLayer, const Real implicit_fac, const bool use_mrf_countergradient)

template<int stagdir>
void	ImplicitDiffForMomLU_T (const Box &bx, const Box &domain, const int level, const double dt_d, const Array4< const int > &col_kext, const Array4< const Real > &cell_data, const Array4< Real > &face_data, const Array4< const Real > &tau, const Array4< const Real > &tau_corr, const Array4< const Real > &z_nd, const Array4< const Real > &detJ, const GpuArray< Real, AMREX_SPACEDIM > &cellSizeInv, const Array4< const Real > &mu_turb, const SolverChoice &solverChoice, const BCRec *bc_ptr, const bool use_SurfLayer, const Real implicit_fac, const bool use_ysu_mom_countergradient)

Macro Definition Documentation

◆ INSTANTIATE_IMPLICIT_DIFF_FOR_MOM_LU

#define INSTANTIATE_IMPLICIT_DIFF_FOR_MOM_LU ( STAGDIR )

Value:

    template void ImplicitDiffForMomLU_T<STAGDIR> ( \
        const Box&, \
        const Box&, \
        const int, \
        const double, \
        const Array4<const int >&, \
        const Array4<const Real>&, \
        const Array4<      Real>&, \
        const Array4<const Real>&, \
        const Array4<const Real>&, \
        const Array4<const Real>&, \
        const Array4<const Real>&, \
        const GpuArray<Real, AMREX_SPACEDIM>&, \
        const Array4<const Real>&, \
        const SolverChoice&, \
        const BCRec*, \
        const bool, \
        const Real, \
        const bool);

Function Documentation

◆ ImplicitDiffForMomLU_T()

template<int stagdir>

void ImplicitDiffForMomLU_T	(	const Box &	bx,
		const Box &	domain,
		const int	level,
		const double	dt_d,
		const Array4< const int > &	col_kext,
		const Array4< const Real > &	cell_data,
		const Array4< Real > &	face_data,
		const Array4< const Real > &	tau,
		const Array4< const Real > &	tau_corr,
		const Array4< const Real > &	z_nd,
		const Array4< const Real > &	detJ,
		const GpuArray< Real, AMREX_SPACEDIM > &	cellSizeInv,
		const Array4< const Real > &	mu_turb,
		const SolverChoice &	solverChoice,
		const BCRec *	bc_ptr,
		const bool	use_SurfLayer,
		const Real	implicit_fac,
		const bool	use_ysu_mom_countergradient
	)

Function for computing the implicit contribution to the vertical diffusion of momentum, over terrain.

This function (explicitly instantiated below) handles staggering in x, y, or z through the template parameter, stagdir. NOTE: implicit diffusion of w has remains an experimental feature and has not been tested yet with terrain.

Parameters

[in]	bx	cell-centered box to loop over
[in]	level	AMR level
[in]	dt_d	time step
[in]	cell_data	conserved cell-centered rho
[in,out]	face_data	conserved momentum
[in]	tau	stress contribution to momentum
[in]	tau_corr	stress contribution to momentum that will be corrected by the implicit solve
[in]	z_nd	nodal array of z
[in]	detJ	Jacobian determinant
[in]	cellSizeInv	inverse cell size array
[in]	mu_turb	turbulent viscosity
[in]	solverChoice	container of parameters
[in]	bc_ptr	container with boundary conditions
[in]	use_SurfLayer	whether we have turned on subgrid diffusion
[in]	implicit_fac	if 1 then fully implicit; if 0 then fully explicit
[in]	use_ysu_mom_countergradient	whether to include YSU momentum countergradient correction

 {
     BL_PROFILE_VAR("ImplicitDiffForMom_T()",ImplicitDiffForMom_T);
  
     Real dt = static_cast<Real>(dt_d);
  
     // setup quantities for getRhoAlphaAtFaces()
     DiffChoice dc = solverChoice.diffChoice;
     TurbChoice tc = solverChoice.turbChoice[level];
     bool l_consA  = (dc.molec_diff_type == MolecDiffType::ConstantAlpha);
     bool l_turb   = tc.use_kturb;
     // The off-diagonal correction strains for u/v contain a factor of 1/2,
     // while the diagonal correction strain for w does not.
     constexpr Real molec_fac = (stagdir == 2) ? two : one;
     Real mu_eff = (l_consA) ? molec_fac * dc.dynamic_viscosity / dc.rho0_trans
                             : molec_fac * dc.dynamic_viscosity;
  
     // g(S*) coefficient
     // stagdir==0: tau_corr = myhalf * du/dz * mu_tot
     // stagdir==1: tau_corr = myhalf * dv/dz * mu_tot
     // stagdir==2: tau_corr =       dw/dz * mu_tot
     constexpr Real gfac = (stagdir == 2) ? two/three : one;
  
     // offsets used to average to faces
     constexpr int ioff = (stagdir == 0) ? 1 : 0;
     constexpr int joff = (stagdir == 1) ? 1 : 0;
  
     // Box bounds
     int ilo = bx.smallEnd(0);
     int ihi = bx.bigEnd(0);
     int jlo = bx.smallEnd(1);
     int jhi = bx.bigEnd(1);
     int box_klo = bx.smallEnd(2);
     int box_khi = bx.bigEnd(2);
     amrex::ignore_unused(ilo, ihi, jlo, jhi);
  
     // Vertical staggering of this component, and the vertical extent of the domain for the
     // faces we are solving on.  A column whose top or bottom is not the domain boundary ends
     // at a coarse/fine boundary instead, and must not be given a physical BC row.
     constexpr int knodal = (stagdir == 2) ? 1 : 0;
     const int dom_klo = domain.smallEnd(2);
     const int dom_khi = domain.bigEnd(2) + knodal;
  
     // Temporary FABs for tridiagonal solve (allocated on column)
     //   A[k] * x[k-1] + B[k] * x[k] + C[k+1] = RHS[k]
     amrex::FArrayBox RHS_fab, soln_fab, coeffG_fab;
            RHS_fab.resize(bx,1, amrex::The_Async_Arena());
           soln_fab.resize(bx,1, amrex::The_Async_Arena());
         coeffG_fab.resize(bx,1, amrex::The_Async_Arena());
     auto const& RHS_a        =        RHS_fab.array();
     auto const& soln_a       =       soln_fab.array();
     auto const& coeffG_a     =     coeffG_fab.array();
  
     Real dz_inv = cellSizeInv[2];
  
     int bc_comp = BCVars::xvel_bc + stagdir;
     bool ext_dir_on_zlo  = (bc_ptr[bc_comp].lo(2) == ERFBCType::ext_dir ||
                             bc_ptr[bc_comp].lo(2) == ERFBCType::ext_dir_prim);
     bool ext_dir_on_zhi  = (bc_ptr[bc_comp].hi(2) == ERFBCType::ext_dir ||
                             bc_ptr[bc_comp].hi(2) == ERFBCType::ext_dir_prim);
     bool foextrap_on_zlo = (bc_ptr[bc_comp].lo(2) == ERFBCType::foextrap);
     bool foextrap_on_zhi = (bc_ptr[bc_comp].hi(2) == ERFBCType::foextrap);
     amrex::ignore_unused(foextrap_on_zlo,foextrap_on_zhi);
  
     AMREX_ASSERT_WITH_MESSAGE(foextrap_on_zlo || ext_dir_on_zlo || use_SurfLayer,
                               "Unexpected lower BC for momentum used with implicit vertical diffusion");
     AMREX_ASSERT_WITH_MESSAGE(foextrap_on_zhi || ext_dir_on_zhi,
                               "Unexpected upper BC for momentum used with implicit vertical diffusion");
  
     Real Fact = implicit_fac * dt * dz_inv;
  
 #ifdef AMREX_USE_GPU
     ParallelFor(makeSlab(bx,2,0), [=] AMREX_GPU_DEVICE (int i, int j, int)
     {
 #else
     for (int j(jlo); j<=jhi; ++j) {
       for (int i(ilo); i<=ihi; ++i) {
 #endif
  
           // The face column at (i,j) is shared by the two cell columns on either side of it,
           // and so by the box on either side of a grid seam.  Solve only over the vertical
           // range that both of those cell columns cover: the boxes on the two sides then
           // build the same tridiagonal system and get the same answer, so the duplicated
           // faces along the seam stay consistent.  Faces outside that range sit on a lateral
           // coarse/fine boundary and are set from the coarse level, so it is consistent to
           // leave them to the explicit update.
           const int klo = amrex::max(box_klo, col_kext(i-ioff,j-joff,0,0),
                                               col_kext(i     ,j     ,0,0));
           const int khi = amrex::min(box_khi, col_kext(i-ioff,j-joff,0,1)+knodal,
                                               col_kext(i     ,j     ,0,1)+knodal);
  
           // Fewer than two faces in common -- nothing to invert, leave the column explicit
 #ifdef AMREX_USE_GPU
           if (khi <= klo) { return; }
 #else
           if (khi <= klo) { continue; }
 #endif
           // Notes:
           //
           // - In DiffusionSrcForMom (e.g., for x-mom)
           //
           //     Real diffContrib = ...
           //                      + (tau13(i,j,k+1) - tau13(i,j,k)) / dzinv
           //     rho_u_rhs(i,j,k) -= diffContrib;  // note the negative sign
           //
           // - We need to scale the explicit _part_ of `tau13` (for x-mom) by (1 - implicit_fac)
           //   The part that needs to be scaled is stored in `tau_corr`.
           //   E.g., tau13 = myhalf * (du/dz + dw/dx)
           //         tau13_corr = myhalf * du/dz
           //
           // - The momentum (`face_data`) was set to `S_old + S_rhs * dt`
           //   prior to including "ERF_Implicit.H". Recall that S_rhs includes
           //   sources from advection and other forcings, not just diffusion.
           //
           // - To correct momentum, we need to subtract `implicit_fac * diffContrib_corr`
           //   from S_rhs to recover `(1 - implicit_fac) * diffContrib_corr`,
           //   where `diffContrib_corr = -d(tau_corr)/dz`. The negative sign
           //   comes from our convention for the RHS diffusion source.
           //
           //   Subtracting a negative gives the += below; multiply by dt to
           //   get the intermediate momentum on the RHS of the tridiagonal
           //   system.
           //
           // - With a surface_layer BC, tau13/23 holds the vertical flux -d_z(k*u_i)
           //   directly. We must use tau at klo (not tau_corr) with SL BCs.
           //
           // - The detJ for divergence was multiplied through.
           //   Therefore, it doesn't show up in the A/B denominator,
           //   but it does modify B and the RHS.
           //
           // - Finally, the terms ~ RHS += (tau_corr_hi - tau_corr_lo) / dz (below)
           //   essentially undo the explicit diffusion update that will be
           //   handled here implicitly.
  
           // Bottom boundary coefficients and RHS for L decomp
           //===================================================
           Real rhoface, rhoAlpha_lo, rhoAlpha_hi;
           Real detJface, met_h_zeta_lo, met_h_zeta_hi;
           Real a_tmp, b_tmp, c_tmp, inv_b2_tmp;
           {
               detJface = myhalf * (detJ(i,j,klo) + detJ(i-ioff,j-joff,klo));
               rhoface  = myhalf * (cell_data(i,j,klo,Rho_comp) + cell_data(i-ioff,j-joff,klo,Rho_comp));
               getRhoAlphaForFaces(i, j, klo, ioff, joff, rhoAlpha_lo, rhoAlpha_hi,
                                   cell_data, mu_turb, mu_eff,
                                   l_consA, l_turb);
  
               met_h_zeta_lo = myhalf * ( Compute_h_zeta_AtKface(i     ,j     ,klo  ,cellSizeInv,z_nd)
                                        + Compute_h_zeta_AtKface(i-ioff,j-joff,klo  ,cellSizeInv,z_nd) );
               met_h_zeta_hi = myhalf * ( Compute_h_zeta_AtKface(i     ,j     ,klo+1,cellSizeInv,z_nd)
                                        + Compute_h_zeta_AtKface(i-ioff,j-joff,klo+1,cellSizeInv,z_nd) );
  
               a_tmp = zero;
               c_tmp = -Fact * gfac * rhoAlpha_hi * dz_inv / met_h_zeta_hi;
  
               RHS_a(i,j,klo) = detJface * face_data(i,j,klo); // NOTE: this is momenta; solution is velocity
  
               // BCs: Dirichlet (u_i = val), slip wall (w = 0), or surface layer (w = 0) --
               //      but only where the bottom of this column really is the domain boundary
               if (klo != dom_klo) {
                   // Coarse/fine (or grid) boundary: close the system with the known face value
                   // below instead of imposing a physical BC at an interior height
                   a_tmp = -Fact * rhoAlpha_lo * dz_inv / met_h_zeta_lo;
                   RHS_a(i,j,klo) += Fact * gfac * (tau_corr(i,j,klo+1) - tau_corr(i,j,klo));
                   const Real rho_below = myhalf * ( cell_data(i     ,j     ,klo-1,Rho_comp)
                                                   + cell_data(i-ioff,j-joff,klo-1,Rho_comp) );
                   RHS_a(i,j,klo) -= a_tmp * (face_data(i,j,klo-1) / rho_below);
               } else if (ext_dir_on_zlo) {
                   RHS_a(i,j,klo) += Fact * gfac * (tau_corr(i,j,klo+1) - tau_corr(i,j,klo));
                   if (stagdir==2) {
                       c_tmp = zero;
                       RHS_a(i,j,klo) = zero;
                   } else {
                       // NOTE: wall is 1/2 dz away (2 dz_inv)
                       a_tmp = -two * Fact * rhoAlpha_lo * dz_inv / met_h_zeta_lo;
                       const Real rho_wall = myhalf * ( cell_data(i     ,j     ,klo-1,Rho_comp)
                                                      + cell_data(i-ioff,j-joff,klo-1,Rho_comp) );
                       const Real wall_velocity = face_data(i,j,klo-1) / rho_wall;
                       RHS_a(i,j,klo) -= a_tmp * wall_velocity;
                   }
               } else if (use_SurfLayer) {
                   // NOTE: tau = -mu*d_z(u_i) w/ SL
                   RHS_a(i,j,klo) += Fact * gfac * (tau_corr(i,j,klo+1) - tau(i,j,klo));
                   RHS_a(i,j,klo) += Fact * tau(i,j,klo);
               } else {
                   // NOTE: FOEXTRAP has zero lower flux (nothing to add to RHS)
                   RHS_a(i,j,klo) += Fact * gfac * (tau_corr(i,j,klo+1) - tau_corr(i,j,klo));
               }
  
               // Add YSU momentum countergradient correction at bottom boundary.
               // NOTE: As for the scalars, the lower face at klo carries no
               //       countergradient flux -- the surface stress is supplied by the
               //       surface layer model or the wall BC above.  Only the upper face
               //       contributes here.
               // NOTE: The sign matches the scalar path: tau_i3 = -rho*K*(du_i/dz - gamma_i),
               //       so the countergradient piece of the flux is +rho*K*gamma_i and its
               //       divergence enters the RHS with a minus sign.
               if (use_ysu_mom_countergradient && stagdir < 2) {
                   const int hgam_comp = (stagdir == 0) ? EddyDiff::HGAMU_v : EddyDiff::HGAMV_v;
                   // Average HGAM* to the staggered face
                   const Real gam_klo = myhalf * (mu_turb(i,j,klo  ,hgam_comp) + mu_turb(i-ioff,j-joff,klo  ,hgam_comp));
                   const Real gam_kp1 = myhalf * (mu_turb(i,j,klo+1,hgam_comp) + mu_turb(i-ioff,j-joff,klo+1,hgam_comp));
                   const Real gam_hi  = myhalf * (gam_klo + gam_kp1);
                   RHS_a(i,j,klo) -= Fact * rhoAlpha_hi * gam_hi;
               }
  
               b_tmp      = detJface * rhoface - a_tmp - c_tmp;
               inv_b2_tmp = one;
  
               RHS_a(i,j,klo)    /= b_tmp;         // NOTE: this is now "rho"
               coeffG_a(i,j,klo)  = c_tmp / b_tmp; // NOTE: this is now "gamma"
           }
  
           // Build the coefficients and RHS for L decomp
           //===================================================
           for (int k(klo+1); k < khi; k++) {
               detJface = myhalf * (detJ(i,j,k) + detJ(i-ioff,j-joff,k));
               rhoface  = myhalf * (cell_data(i,j,k,Rho_comp) + cell_data(i-ioff,j-joff,k,Rho_comp));
               getRhoAlphaForFaces(i, j, k, ioff, joff, rhoAlpha_lo, rhoAlpha_hi,
                                   cell_data, mu_turb, mu_eff,
                                   l_consA, l_turb);
  
               met_h_zeta_lo = myhalf * ( Compute_h_zeta_AtKface(i     ,j     ,k  ,cellSizeInv,z_nd)
                                        + Compute_h_zeta_AtKface(i-ioff,j-joff,k  ,cellSizeInv,z_nd) );
               met_h_zeta_hi = myhalf * ( Compute_h_zeta_AtKface(i     ,j     ,k+1,cellSizeInv,z_nd)
                                        + Compute_h_zeta_AtKface(i-ioff,j-joff,k+1,cellSizeInv,z_nd) );
  
               a_tmp      = -Fact * rhoAlpha_lo * dz_inv / met_h_zeta_lo;
               c_tmp      = -Fact * rhoAlpha_hi * dz_inv / met_h_zeta_hi;
               b_tmp      = detJface * rhoface - a_tmp - c_tmp;
               inv_b2_tmp = one / (b_tmp - a_tmp * coeffG_a(i,j,k-1));
  
               RHS_a(i,j,k)    = detJface * face_data(i,j,k); // NOTE: this is momenta; solution is velocity
               RHS_a(i,j,k)   += Fact * gfac * (tau_corr(i,j,k+1) - tau_corr(i,j,k));
  
               // Add YSU momentum countergradient correction
               if (use_ysu_mom_countergradient && stagdir < 2) {
                   const int hgam_comp = (stagdir == 0) ? EddyDiff::HGAMU_v : EddyDiff::HGAMV_v;
                   // Average HGAM* to the staggered face
                   const Real gam_k   = myhalf * (mu_turb(i,j,k  ,hgam_comp) + mu_turb(i-ioff,j-joff,k  ,hgam_comp));
                   const Real gam_km1 = myhalf * (mu_turb(i,j,k-1,hgam_comp) + mu_turb(i-ioff,j-joff,k-1,hgam_comp));
                   const Real gam_kp1 = myhalf * (mu_turb(i,j,k+1,hgam_comp) + mu_turb(i-ioff,j-joff,k+1,hgam_comp));
                   const Real gam_hi  = myhalf * (gam_k + gam_kp1); // at k+1/2
                   const Real gam_lo  = myhalf * (gam_k + gam_km1); // at k-1/2
                   RHS_a(i,j,k) -= Fact * (rhoAlpha_hi * gam_hi - rhoAlpha_lo * gam_lo);
               }
  
               RHS_a(i,j,k)    = (RHS_a(i,j,k) - a_tmp * RHS_a(i,j,k-1)) * inv_b2_tmp; // NOTE: This is now "rho"
               coeffG_a(i,j,k) = c_tmp * inv_b2_tmp; // NOTE: this is now "gamma"
           } // k
  
           // Top boundary coefficients and RHS for L decomp
           //===================================================
           {
               detJface = myhalf * (detJ(i,j,khi) + detJ(i-ioff,j-joff,khi));
               rhoface  = myhalf * (cell_data(i,j,khi,Rho_comp) + cell_data(i-ioff,j-joff,khi,Rho_comp));
               getRhoAlphaForFaces(i, j, khi, ioff, joff, rhoAlpha_lo, rhoAlpha_hi,
                                   cell_data, mu_turb, mu_eff,
                                   l_consA, l_turb);
  
               met_h_zeta_lo = myhalf * ( Compute_h_zeta_AtKface(i     ,j     ,khi  ,cellSizeInv,z_nd)
                                        + Compute_h_zeta_AtKface(i-ioff,j-joff,khi  ,cellSizeInv,z_nd) );
               met_h_zeta_hi = myhalf * ( Compute_h_zeta_AtKface(i     ,j     ,khi+1,cellSizeInv,z_nd)
                                        + Compute_h_zeta_AtKface(i-ioff,j-joff,khi+1,cellSizeInv,z_nd) );
  
               a_tmp = -Fact * gfac * rhoAlpha_lo * dz_inv / met_h_zeta_lo;
               c_tmp = zero;
  
               RHS_a(i,j,khi)  = detJface * face_data(i,j,khi); // NOTE: this is momenta; solution is velocity
               RHS_a(i,j,khi) += Fact * gfac * (tau_corr(i,j,khi+1) - tau_corr(i,j,khi));
  
               // BCs: Dirichlet (u_i = val), slip wall (w = 0) --
               //      but only where the top of this column really is the domain boundary
               if (khi != dom_khi) {
                   // Coarse/fine (or grid) boundary: close the system with the known face value
                   // above instead of imposing a physical BC at an interior height
                   c_tmp = -Fact * rhoAlpha_hi * dz_inv / met_h_zeta_hi;
                   const Real rho_above = myhalf * ( cell_data(i     ,j     ,khi+1,Rho_comp)
                                                   + cell_data(i-ioff,j-joff,khi+1,Rho_comp) );
                   RHS_a(i,j,khi) -= c_tmp * (face_data(i,j,khi+1) / rho_above);
               } else if (ext_dir_on_zhi) {
                   if (stagdir==2) {
                       a_tmp = zero;
                       RHS_a(i,j,khi) = zero;
                   } else {
                       // NOTE: wall is 1/2 dz away (2 dz_inv)
                       c_tmp = -two * Fact * rhoAlpha_hi * dz_inv / met_h_zeta_hi;
                       const Real rho_wall = myhalf * ( cell_data(i     ,j     ,khi+1,Rho_comp)
                                                      + cell_data(i-ioff,j-joff,khi+1,Rho_comp) );
                       const Real wall_velocity = face_data(i,j,khi+1) / rho_wall;
                       RHS_a(i,j,khi) -= c_tmp * wall_velocity;
                   }
               }
  
               b_tmp      = detJface * rhoface - a_tmp - c_tmp;
               inv_b2_tmp = one / (b_tmp - a_tmp * coeffG_a(i,j,khi-1));
  
               // First solve
               soln_a(i,j,khi) = (RHS_a(i,j,khi) - a_tmp * RHS_a(i,j,khi-1)) * inv_b2_tmp;
           }
  
           // Back sweep the U decomp solution
           //===================================================
           for (int k(khi-1); k>=klo; --k) {
               soln_a(i,j,k) = RHS_a(i,j,k) - coeffG_a(i,j,k) * soln_a(i,j,k+1);
           }
  
           // Convert back to momenta
           //===================================================
           for (int k(klo); k<=khi; ++k) {
               rhoface = myhalf * (cell_data(i,j,k,Rho_comp) + cell_data(i-ioff,j-joff,k,Rho_comp));
               face_data(i,j,k) = rhoface * soln_a(i,j,k);
           }
  
 #ifdef AMREX_USE_GPU
     });
 #else
       } // i
     } // j
 #endif
 }

Here is the call graph for this function:

◆ ImplicitDiffForStateLU_T()

void ImplicitDiffForStateLU_T	(	const Box &	bx,
		const Box &	domain,
		const int	level,
		const int	n,
		const double	dt_d,
		const GpuArray< Real, AMREX_SPACEDIM *2 > &	bc_neumann_vals,
		const Array4< Real > &	cell_data,
		const Array4< const Real > &	z_nd,
		const Array4< const Real > &	detJ,
		const GpuArray< Real, AMREX_SPACEDIM > &	cellSizeInv,
		const Array4< const Real > &	scalar_zflux,
		const Array4< const Real > &	mu_turb,
		const SolverChoice &	solverChoice,
		const BCRec *	bc_ptr,
		const bool	use_SurfLayer,
		const Real	implicit_fac,
		const bool	use_mrf_countergradient
	)

Function for computing the implicit contribution to the vertical diffusion of theta, with terrain.