sandbox/Antoonvh/diurnalSCM.c
Note: This page is best viewed with the Firefox web browser.
A Single-Column model for the Atmospheric Diurnal Cycle
We time integrate a 1D evolution equation for the vertical profiles of the horizontal velocity components (u,v) and the buoyancy (b), using an adaptive bitree grid, a generic timeloop iterator and the reaction-diffusion solver.
#include "grid/bitree.h"
#include "run.h"
#include "diffusion.h"
For the turbulent mixing we use the following mixing functions (f(Ri)), for stable (fris
) and unstable (friu
) stratifications,
//#define fris(Ri) (sq((1 - (Ri/0.20)))*(Ri < 0.20)) // Critical Ri
#define fris(x) (exp(-10.*x)) // Exponential
//#define fris(x) (1./(1. + (10*x*(1.+8.*x)))) // Long Tail
#define friu(Ri) (sqrt(1. - (18.*Ri))) // Holtslag en Boville 1992
// Note: no surface f(Ri)
The thermodynamic variable is forced at the surface with the surface buoyancy flux (B) according to a a simple energy balance:
\displaystyle B = Q* + G,
where Q* in a prescribed netto radiation and G is a ‘feed back flux’ that is dynamically evaluated as:
\displaystyle G = -\Lambda \left( b_{surf} - b_{ref} \right),
with b_{ref} a constant reference temperature/buoyancy scale and \Lambda a coupling strength. The buoyancy at the surface (b_{surf}) is evaluated using the lowest two cell values for b (b_1, b_2) and a sub-grid-scale model that assumes a logaritmic profile for b(z).
\displaystyle b_{surf} = \frac{b_2 - b_1*c}{1-c},
with,
\displaystyle c = \frac{\mathrm{ln}\left(\frac{4\Delta}{z_{0,h}}\right)}{\mathrm{ln}\left(\frac{\Delta}{z_{0,h}}\right)}
that is the result from an integration exercise, assuming that the grid-cell size (\Delta) is much larger than the roughness length for heat (z_{0,h}); \Delta \gg z_{0,h}. The surface buoyancy can be evaluated from the lowest to cells according to a lookup table that stores the values for c(\Delta). Furthermore, we #define
macros to compute G(GFLX
) and Q*(Qn
), such that we can readily evaluate B.
#define BSURF ((b[1] - b[]*c[level])/(1. - c[level]))
#define GFLX (-Lambda*(BSURF - bo))
Prescribing,
\displaystyle Q*= \mathrm{max}\left[ B_1\mathrm{sin}\left(\frac{2 \pi t}{T} \right), B_1 \right].
#define Qn (max(B0*sin(2.*M_PI*t/T), B1))
We declare fields for u, v and b as u
, v
and b
, respectively and also set the value for L_c = 1030m and T = 24\times3600s.
scalar u[], v[], b[];
double Lc = 1030.;
double T = 24.*3600.;
global parameters are declared for \Lambda, B_0, B_1, N, U_{geo}, z_{0,m} and z_{0,h}.
double Lambda, B0, B1, f, Nv, Ugeo, zom, zoh;
The numerical values of which will be computed from the values of the dimensionless groups:
double Pi1 = -6;
double Pi2 = 2160.;
double Pi3 = 10.;
double Pi4 = 5366.;
double Pi5; // This one will be varied
double Pi6 = 5150.;
double Pi7 = 1.;
We initialize/declare some usefull variables for the deep soil
buoyancy, the Von Karmann constant, a lookup table to help determine b_{surf} and the surface friction velocity (u*), the Mixed-layer depth (\mathcal{L}_c), the number of averaging iterations (for the inversion data), the maximum level of refinement, a structure for the multigrid-solver details, the maximum windspeed, the inversion, U at \approx40m and a file pointer for the output of the inversion strengths.
double bo = 0., k = 0.4;
double c[20], lut[20], Lmix;
int ni, maxlevel = 9;
mgstats mgb;
double Umax, inv, u40;
FILE * fpi;
We set boundary conditions where left
and right
are to be interpreted as the surface and the domain top, respectively.
b[left] = dirichlet (BSURF);
b[right] = dirichlet (bo + sq(Nv)*x);
u[left] = dirichlet (0.);
v[left] = dirichlet (0.);
int main(){
The file is opened and the parameter values are calculated from the dimensionless groups.
fpi = fopen ("Inversions", "w");
Nv = Pi2/T; //Initial stratification strength
f = Pi3/T; //Coriolis parameter
B0 = sq(Lc*Nv)*M_PI/(2.*T); //B_0
B1 = B0/Pi1; //B_1
Lambda = sqrt(B0*T)/Pi4; //Lambda
zom = Lc / Pi6; //z_0,m
zoh = zom * Pi7; //z_0,h
L0 = 3.*Lc; //Domain size
We vary the value of \Pi_5 from 1 to 8 in 57 runs and compute the value of U_{geo} accordingly.
for (Pi5 = 1.; Pi5 <= 8; Pi5 += 0.125){
Ugeo = Pi5*pow(B0*Lc, 1./3.);
init_grid (128);
run();
}
}
This event initializes the setup.
First, the inversion-strength data is reset.
ni = 0;
inv = Umax = u40 = 0.;
For accurate time integration we initially set a small timestep (1 second) and only allow a small tolerance on the residual for the Poisson problem. The timestep will be adaptive, the TOLERANCE
remains small.
DT = T/(3600.*24.);
TOLERANCE = 1E-4;
Solution fields are initialized after the near-surface grid is refined to the maximum resolution.
refine(x < (Lc/10.) && level < maxlevel);
foreach(){
u[] = Ugeo;
v[] = 0.;
b[] = bo + sq(Nv)*x;
}
For the computations of b_{surf}(c
) and the friction velocity u* (lut
), we create lookup tables for the various possible grid resolutions thay may be used at the surface.
for (int j = 1; j <= maxlevel; j++){
double d = (L0/((double)(1 << j)))/zoh;
double d2 = (L0/((double)(1 << j)))/zom;
c[j] = (log(4.*d) - 1.)/(log(d) - 1.);
lut[j] = sq(k/(log(d2) - 1.));
}
}
Time integration
We time integrate the evolution equation of the atmospheric profiles:
We allocate temporary fields to store the tendency terms and the diffusivity.
scalar rx[],ry[],rb[];
face vector kh[];
double B = 0;
double ws = 0;
The momentum is forced with a height-and-time constant horizontal pressure gradient (-\nabla P) and is affected by back ground rotation according to the Coriolis Parameter (f). Introducing,
\displaystyle \overrightarrow{U_{geo}} = \frac{\overrightarrow{k}}{\rho f}\times \nabla P,
that is known as the geostrophic wind.
foreach(){
rx[] = f*v[];
ry[] = f*(Ugeo - u[]);
rb[] = 0.;
In the bottom cell, we add the tendency due to the buoyancy flux.
if (x < Delta){
B = (Qn + GFLX);
rb[] += B / Delta;
For the momentum flux at the surface, we use a robust law-off-the-wall in the lowest grid cell (z+z_0 < \Delta).
\displaystyle u(z) = \frac{u*}{\kappa}\mathrm{ln}\left(\frac{z}{z_{0.m}}\right)
with u* the friction velocity, \kappa=0.4 the Von Karmann constant, \mathrm{ln}(x) is natural logarithm of a dummy variable x and z_{0,m} the roughness length for momentum. We can express a tendency in the lowest cell [u] due to the surface friction:
\displaystyle u*^2 = \frac{u}{\|u\|} \left( \frac{uk}{\mathrm{ln} \left(\frac{\Delta}{z_{0,m}} \right)} \right)^2.
The corresponding prefactors are already computed and stored in the lut
array.
rx[] += -sign(u[])*lut[level]*sq(u[])/Delta;
ry[] += -sign(v[])*lut[level]*sq(v[])/Delta;
}
}
The local eddy diffusivity (K) is written as,
\displaystyle K = lV_*
where l is the mixing length,
\displaystyle l = \mathrm{min} \left[ kz, 70m \right],
with k= 0.4 the Von Karman constant. V_* is the mixing velocity scale due to shear and convection, inspired by Troen en Mahrt (1986) we write,
\displaystyle V_* = \sqrt{w_c^2 + \left(lSf(Ri)\right)^2}
with w_c the vertical velocity variance for convection (B>0 and z < \mathcal{L}_c),
\displaystyle w_c =3 w_d z/h \left( 1-z/\mathcal{L}_c \right) ^ 2,
with w_d the Deardorf velocity scale,
\displaystyle w_d = (B\mathcal{L}_c)^{1/3},
and \mathcal{L}_c the height of the well mixed layer,
\displaystyle \mathcal{L}_c =\sqrt{\frac{2\int b - N^2z \mathrm{d}z}{N^2}}
for which we followed the works of Van Heerwaarden, Mellado etc.
Noting that w_c = 0 for the times with Q*(t) < 0 and heights z > \mathrm{L}_c.
The term lSf(Ri) is computed conforming to the definitions in Van Hooft et al. (2018b).
Lmix = 0.;
foreach()
Lmix += (b[] - x*sq(Nv)) * Delta;
if (Lmix > 0.)
Lmix = sqrt (Lmix*2./sq(Nv));
if (B > 0 && Lmix > 0)
ws = pow(B*Lmix, 1./3.);
foreach_face(){
double sqd = (sq((u[] - u[-1])/(Delta)) + sq((v[] - v[-1])/(Delta)));
double Ri = ((b[] - b[-1])/(Delta))/(sqd + 0.00001);
double fRi;
if (Ri < 0)
fRi = friu(Ri);
else
fRi = fris(Ri);
double l = min(k*x, (70./1030.)*Lc); // The mixing length
double fraction = 0;
if (Qn > 0)
fraction = 3.*(x/Lmix * sq(1. - x/Lmix))*(x < Lmix);
double Vs = sqrt(fraction*sq(ws) + sq(l)*sqd*sq(fRi));
kh.x[] = l*Vs;
}
boundary(all);
dt = dtnext(DT);
Now we have all the ingredients for the reaction-diffusion problem.
int n = 0;
mgb = diffusion (u, dt, kh, r = rx);
n += mgb.i;
mgb = diffusion (v, dt, kh, r = ry);
n += mgb.i;
mgb = diffusion (b, dt, kh, r = rb);
n += mgb.i;
Based on the convergence properties we adapt the timestep.
if (n > 10) //Quickly reduce the timestep if things get rough
DT = max (DT/(1+((double)n/10.)), T/(24.*3600.));
if (n < 5) //Slowly increase the timestep when time integration is easy.
DT = min (DT*(1+((double)n/100.)), 20.*T/(24.*3600.));
}
The grid is adapted based on the wavelet-estimated error for the discretized representation of the solution field. The criteria are chosen to be \zeta_{u,v} = U_{geo}/20 and \zeta_b = b_{c,\Lambda}/50. These values work well enough.
event adapt (i++){
double ue = Ugeo/20.;
double be = sq(Nv)*Lc/50.;
adapt_wavelet({u,v,b}, (double[]){ue, ue, be}, maxlevel);
}
Output
Output consists of instantanious porofiles that are outputted 15 times per phyisical hour.
event profile(t += T/(24*15)){
char fname[99];
sprintf (fname, "ProfsPi5%g", Pi5);
static FILE * fp = fopen (fname, "w");
foreach(){
if (x < Delta)
fprintf (fp, "0 0 0 %g \n", BSURF);
fprintf (fp, "%g %g %g %g \n", x, u[], v[], b[]);
}
}
We also monitor the occurance of a low level jet. The numerical values are not very robust for the various mixing function (f(\mathrm{Ri})) formulations, and hence not presented in the associated work.
double xuvb[4][1000] = {0};
double trec = 0;
event llj (t = T/2; i += 5){
bool new_record = false;
foreach()
if (sq(u[]) + sq(v[]) > sq(Umax))
new_record = true;
if (new_record){
trec = t;
int j = 0;
foreach(){
xuvb[0][j] = x;
xuvb[1][j] = u[];
xuvb[2][j] = v[];
xuvb[3][j++] = b[];
}
xuvb[0][j] = -1.;
}
}
Each timestep, we monitor the buoyancy difference (i.e. the inversion) between z = 0 and z = L_c/20 and the windspeed at that height. These are summed up over the 19th hour of simulation.
event inversion (t = T*18./24.; t<= T*19./24.; i++){
double zp1 = 0.;
double zp2 = Lc / 20;
inv += interpolate (b, zp2) - interpolate (b, zp1);
u40 += sqrt(sq(interpolate (u, zp2)) + sq(interpolate (v, zp2)));
ni++;
}
Each ten timesteps we output some interesting solution diagnostics so that we can see how the solution has evolved over time. We output (Q*(t), G(t), B(t), b_{surf}, b_{mix} and \mathcal{L}_c to a file named: timeseriesPi5
followed by a \Pi_5-value identifier.
event timeseries (i += 10){
char fname[99];
sprintf (fname, "timeseriesPi5%g", Pi5);
static FILE * fp = fopen (fname, "w");
double B, G, bs, bmix;
Lmix = 0;
foreach()
Lmix += (b[] - x*sq(Nv)) * Delta;
if (Lmix > 0.)
Lmix = sqrt (Lmix*2./sq(Nv));
bmix = Lmix * sq(Nv);
foreach_boundary (left){
G = GFLX;
B = Qn + G;
bs = BSURF;
}
if (i == 0)
fprintf (fp, "t\tQn\tG\tB\tbs\tbmix\tLmix\n");
fprintf (fp, "%g\t%g\t%g\t%g\t%g\t%g\t%g\n",
t, Qn, G, B, bs, bmix, Lmix);
}
The last event
The run()
stops once t = T and then we output the (iteration) averaged inversion data and the low-level jet data.
event stop (t = T){
int m = 0;
char fname[99];
sprintf (fname, "Profjet%g",Pi5);
static FILE * fpj = fopen (fname, "w");
while (xuvb[0][m] != -1.){
fprintf(fpj, "%g\t%g\t%g\t%g\n",
xuvb[0][m], xuvb[1][m],xuvb[2][m], xuvb[3][m]);
m++;
}
inv /= (double)ni;
u40 /= (double)ni;
fprintf (fpi, "%g\t%g\t%g\n", Pi5, inv, u40);
}
A movie
On systems with gnuplot
and ffmpeg
installed, we can choose to generate animations of the simulation results:
#define GNUPLOT_AND_FFMPEG 1 //Switch for movie output
#if GNUPLOT_AND_FFMPEG
We initialize a pipeline for the plots that will later be turned into a movie.
static FILE * gnuplotPipe;
event init(t = 0){
gnuplotPipe = popen ("gnuplot", "w");
fprintf (gnuplotPipe,
"set term pngcairo size 1200,600 enhanced font 'Times ,18'\n"
"set yr [0: 1200]\n"
"set ylabel 'height [m]'\n"
"set grid\n"
"set size square\n");
}
We output 15 plots per physical hour.
int frame = 0;
event plot (t += T/(24*15)){
if (fabs(Pi5 - round(Pi5)) < 0.001){
fprintf (gnuplotPipe, "set output 'plot%d.png'\n", frame);
fprintf (gnuplotPipe, "set multiplot layout 1,2 title 'Π_5 = %g,"
"time %.02d:%.02d' font 'Times ,25'\n",
Pi5, (int)(t/(3600)), ((int)t%3600)/60);
fprintf (gnuplotPipe,
"set xr [-15: 20]\n"
"set xlabel 'wind speed [m/s]'\n"
"set key top left\n"
);
fprintf (gnuplotPipe, "plot '-' w l lw 5 lt rgb'#11BB11' t 'u',"
"'-' w l lw 5 lt rgb '#BB11BB' t 'v'\n");
foreach()
fprintf (gnuplotPipe, "%g %g\n",u[], x);
fprintf (gnuplotPipe, "e\n");
foreach()
fprintf (gnuplotPipe, "%g %g\n",v[], x);
fprintf (gnuplotPipe, "e\n");
fprintf (gnuplotPipe,
"set xr [260 : 290]\n"
"set key off\n"
"set xlabel 'θ [K]'\n"
);
fprintf (gnuplotPipe, "plot '-' w l lw 5 lt rgb '#CC1111'\n");
Note that the movie displays the potential temperature (\theta) rather than the buoyancy, using \theta_{ref} = 270K and g = 10 ms^{-2}.
foreach()
fprintf (gnuplotPipe, "%g %g\n", b[]*27 + 270, x);
fprintf (gnuplotPipe, "e\n"
"unset multiplot\n");
fflush (gnuplotPipe);
frame++;
}
}
Finally, we render a .mp4
movie from these plots and then remove the plots from the disk.
event moviemaker(t = T){
if (fabs(Pi5 - 8.) < 0.0001){
system ("rm mov_scm.mp4");
system ("ffmpeg -r 25 -f image2 -i plot%d.png -vcodec libx264 mov_scm.mp4");
system ("rm plot*");
}
return 1;
}
#endif