sandbox/Antoonvh/thermo.h

    Moist thermodynamics

    This page implements moist atmospheric thermodynamics, suitable for turbulence-resolving boundarly-layer simulations. In order to simulate clouds in the atmosphere we adopt a model for phase changes. Notice this this happens entirely on the subgrid scale. We take inspirtation from the formulations of DALES and Micro-HH.

    A bunch of global constants are defined. Furthermore, fields for the liquid water potential temperature (\theta_v, thl in Kelvin), the total water specific humidity (q_t, qt in kg/kg) and an approximate thermodynamic pressure (p, pres in Pascal) are declared.

    double Rv = 461.5;        //in J/kg/K
double Rd = 287.;         //in J/kg/K
double L = 2500000.;      //in J/kg
double cpd = 1004.;       //in J/Kg/Kq
#define P00 100000.       //A constant reference pressure
double P0 = 101780., Pe = 7500;//Surface pressure and 1/e length scale
scalar thl[], qt[], pres[];
pres[bottom] = dirichlet (P0);
#ifndef RHO0            
#define RHO0 1.3 //Base density in kg/m^3  
#endif
double rho0 = RHO0;

#if dimension == 1
#define HEIGHT_ABOVE_SURF (x - X0)  
#else
#define HEIGHT_ABOVE_SURF (y - Y0)
#endif

    A sequence of DEFINEs is used to compute the cloud liquid water content as a function of height and temperature. Note that there is no field dedicated to store the liquid water specific humidity.

    #define PRES  (P0*(exp(-HEIGHT_ABOVE_SURF/Pe)))   
#define EXNER (pow(pres[]/P00, Rd/cpd))
#define TK1   (thl[]*EXNER)
#define ES    (610.78*exp(17.27*(TK1 - 273.16)/(TK1 - 35.86)))
#define QSL   (Rd/Rv*(ES/(pres[] - (1. - Rd/Rv)*ES)))
#define QSAT  (QSL*((1. + ((sq(L)/(Rv*cpd*sq(TK1)))*qt[]))/(1 + ((sq(L)/(Rv*cpd*sq(TK1)))*QSL))))
#define QC    (max(qt[] - QSAT, 0))

    The virtual potential temperature in a cell (\theta_v,THV) can be computed from \theta_l, q_c and q_t.

    #define THETA (thl[] + L*QC/(cpd*EXNER)) 
#define THV (THETA*(1. - (1. - Rv/Rd)*qt[] - Rv/Rd*QC))

    The virtual potential temperature is related to the buoyancy,

    double g_acc = 9.81, T_ref = 280.;
#define B_CENTERED (g_acc*(THV - T_ref)/T_ref)

    The thermodynamic pressure

    We need to have a quess at the thermodynamic pressure. It can integrated from the surface pressure (P0) with height according to:

    \displaystyle \frac{\partial \Pi}{\partial z} = -\frac{g}{c_p\theta_{v}} with, \displaystyle \Pi = \left(\frac{p}{p_0} \right)^{\frac{Rd}{c_p}}

    Notice that this pressure definition is not modified by the flow. The assumption being that g is far greater than any acceleration in the flow field. We use a multigrid solver for this purpose.

    #include "integrator.h"
struct STP {
  bool reset;   //Use a guess for the pressure
  double TH;    //The angle (in radians) of the vertical direction
};
int set_thermodynamic_pressure (struct STP stp) {
  if (!stp.TH)
    THETA_ANGLE = -pi/2.; //`Bottom` is the surface
  else
    THETA_ANGLE = stp.TH;
  scalar rhs[], _pres[], Pi[]; 
  double TOL = 0.1;  //Pa
  Pi[bottom] = dirichlet (pow(P0/P00, Rd/cpd));

    Notice that computing the rhs requires the pressure itself. Therefore, we must call the iterative solver iteratively and hope that the pressure converges. As a fist guess we use the current values in pres[], unless the argument reset is true, then we take a typical profile. It is not wise to call this function very often.

      foreach() {
    _pres[] = PRES;
    if (stp.reset)
      pres[]  = _pres[];
    Pi[] = EXNER;
  }
  int j = 0;
  do {
    foreach()
      rhs[] = -g_acc/(cpd*THV);
    integrate_dx (Pi, rhs);
    j++;
    foreach()
      pres[] = pow((Pi[]), cpd/Rd)*P00;
  } while (change (pres, _pres) > TOL && j < NITERMAX);
  if (j == NITERMAX)
    fprintf(ferr, "The thermodynamic pressure was not found within %d cycles\n", j);
  return j;
}

    Alternatively, one may opt for the set_pres () functionality (tip!).

    #include "profile5c.h"
struct SP {
  int C;      // Levels
  bool guess; // Use a guess
};

void set_pres (struct SP sp) {
  int C = (1 << grid->maxdepth);
  if (sp.C)
    C = sp.C;
  double dz = 0.99999*L0/(double)(C - 1); // almost span the entire domain

    First we compute the \theta_v field so we can find its horizontally averaged values.

      scalar thv[];
  foreach() {
    if (sp.guess)
      pres[] = PRES;
    thv[] = THV;
  }
  boundary ({thv});

    A 1D look-up table is created for the pressure (Y), computed via \theta_{v0}, as described in van Heerwaarden et al. (2017).

      double Y[C][3]; //y, thl, p
  memset (Y, 0, sizeof(Y));
  Y[0][0] = Y0 + SEPS;
  Y[0][2] = P0;
  average_over_yp ({thv}, &Y[0][1], Y[0][0]);
  for (int n = 1; n < C; n++) {
    Y[n][0] = Y[n - 1][0] + dz;
    average_over_yp ({thv}, &Y[n][1], Y[n][0]);
    if (pid() == 0)
      Y[n][2] = Y[n - 1][2]*exp(-g_acc*2.*dz/(Rd*pow(Y[n-1][2]/P00, Rd/cpd)*
					      (Y[n][1] + Y[n - 1][1])));
  }

    In parallel, the root broadcasts its data.

    #if _MPI
  MPI_Bcast (&Y[0][0], 3*C, MPI_DOUBLE, 0, MPI_COMM_WORLD);
#endif

    Each cell can find the corresponding pressure via interpolation

      foreach() {
    double H = HEIGHT_ABOVE_SURF;
    int n =  (int)(H/dz) + 1;
    pres[] = ((H - Y[n - 1][0])*Y[n][2] +
	      (Y[n][0] - H)*Y[n - 1][2])/dz;
  }
  boundary ({pres});  
}

    Moist dynamics

    Next, we couple the thermodynamics to the flow. That means q_t and \theta_l become active tracers. It can optionally be swithed off by compiling with a -DNO_DYNAMICS tag.

    #ifndef NO_DYNAMICS
#include "tracer.h"
scalar * tracers = {qt, thl};

event defaults (i = 0) {
  if (is_constant(a.x)) {
    a = new face vector;
    foreach_face()
      a.x[] = 0;
    boundary ((scalar*){a});
  }
  pres.prolongation = refine_linear;
#if TREE
  pres.refine = refine_linear;          //conservative 3rd order vertical interp.
#endif
  thl.gradient = qt.gradient = minmod2; //A slope limiter
}

    The effect of gravity on the flow is described via the buoyancy.

    event acceleration (i++) {
  face vector av = a;
  scalar buoy[];
  foreach()
    if (cm[] > 0)  //?
      buoy[] = B_CENTERED;
  boundary ({buoy});
  foreach_face(y) 
    av.y[] += fm.y[]*(buoy[0,-1] + buoy[])/2.;
}

    In order to find the modified pressures without triggering spurrious currents, a small tolerance is used for a brief while.

    int back_it = 10;
double DEF_TOL = 1e-3, DEF_DT = 60.;

event init (t = 0) {
  DT = min(DT, 1.); 
  TOLERANCE = min(TOLERANCE, 1e-6);
}

event back_to_defaults (i = back_it) {
  TOLERANCE = DEF_TOL;
  DT = DEF_DT; 
}

#endif //NO_DYNAMICS

    extras

    A blue-to-white color bar for cloud visuals

    void cloud (double cmap[NCMAP][3]) {
  for (int i = 0; i < NCMAP; i++) {
    cmap[i][0]= 0.2 + 0.8 * (double)i/(double)NCMAP;
    cmap[i][1]= 0.2 + 0.8 * (double)i/(double)NCMAP;
    cmap[i][2]= 0.9 + 0.1 * (double)i/(double)NCMAP; 
  }
}

    Usage

    Reference

    Heerwaarden, C. C. V., Van Stratum, B. J., Heus, T., Gibbs, J. A., Fedorovich, E., & Mellado, J. P. (2017). MicroHH 1.0: a computational fluid dynamics code for direct numerical simulation and large-eddy simulation of atmospheric boundary layer flows. Geoscientific Model Development, 10(8), 3145-3165.