sandbox/Antoonvh/GABLS2forrev1and2.c
Variations on the second GABLS intercomparison case
This page presents how the mixing length l controls the slope in the observed temperature profile. For a descroption of the set-up see the default set-up
#include "grid/bitree.h"
#include "diffusion.h"
#include "run.h"
#define T1bottom (((((-25*cos(0.22*((t/3600)+16) + 0.2)-10) *(((t/3600)+16)<=17.4))+((((t/3600)+16)<=30)*(((t/3600)+16)>17.4)*((-0.54*((t/3600)+16))+15.2))+((((t/3600)+16)<=41.9)*(((t/3600)+16)>30)*(-7-(25*cos(0.21*((t/3600)+16)+1.8))))+((((t/3600)+16)<=53.3)*(((t/3600)+16)>41.9)*(-0.37*((t/3600)+16)+18.0))+((((t/3600)+16)<=65.6)*(((t/3600)+16)>53.3)*(-4-25*cos(0.22*((t/3600)+16)+2.5)))+((((t/3600)+16)>65.6)*(4.4)))+273.15))
#define maxlevel 9
// Holtslag and Boville with Stable F(Ri) according to vdW
#define fris(Ri) (sq((1-(Ri/0.20)))*(Ri<0.20)) //vdW 2017 Critical Ri
//#define fris(x) (1/(1+(10*x*(1+8*x))))
#define friu(Ri) (sqrt(1-(18.*Ri))) // Holtslag en Boville 1992
#define friubm(Ri,y) ((1-((10*Ri)/(1+75*y*sqrt((x+zo/zo)*fabs(Ri)))))) // Louis 1982
#define friubh(Ri,y) ((1-((15*Ri)/(1+75*y*sqrt((x+zo/zo)*fabs(Ri)))))) // louis 1982
Below some constants are decleared.
double Tref=283.15; // Reference temperture
double Ugeo=-8.; // u-component of the Geostrophic wind
double Vgeo = 3.; // v-component of the Geostrophic wind
double zo=0.01; // Roughness length for heat and momentum
double He= 65.6; //Total hours of simulation
double Rv= 461.5; // Gas constant for water vapor
double Rd = 287.; // Gas constant for dry air
double L= 2500000.; // Latent heat release for evaporization
double cpd= 1004; // Density of air.
double exner;
int nn;
double Up[100];
double Cm,Ch,Cq,qtsat,TK1;
scalar u[],v[],T1[],qt[];
mgstats mgb;
double eu=0.25;
double eb=0.5;
Here the maxim mixing length is declared as a variable lz
so that it can be varied.
double lz = 70;
int j;
int main(){
4 runs.
We run four different runs;
- Default run
- Using a decreased max. mixing length lz
- Using an increased max. mixing length lz
- Default run, without adaptivity.
init_grid(1<<(maxlevel));
L0=4096;
X0=0;
j = 0;
run();
j = 1;
lz=35;
run();
j=2;
lz=140;
run();
j=3;
lz=70;
run();
}
u[left]=dirichlet(0.);
v[left]=dirichlet(0.);
T1[left]=dirichlet(T1bottom);
T1[right]=neumann((20./3000.));
event init(i=0){
exner = pow(972./1000.,Rd/cpd);
TOLERANCE=10E-8;
DT=1;
u.refine=refine_linear;
v.refine=refine_linear;
T1.refine=refine_linear;
foreach(){
u[]=Ugeo;
v[]=Vgeo;
T1[]=(((x<=200))*(288-(2*x/200)))+
((x>200)*(x<=850)*286)+
((x>850)*(x<=900)*(286+2*(x-850)/50))+
((x>900)*(x<=1000)*(288+4*(x-900)/100))+
((x>1000)*(292+20*((x-1000)/3000)));
qt[]=0.0025;
}
if (j!=3){
while(adapt_wavelet({u,v,T1},(double[]){eu,eu,eb},maxlevel,4,{u,v,T1,qt}).nc){
foreach(){
u[]=Ugeo;
v[]=Vgeo;
T1[]=(((x<=200))*(288-(2*x/200)))+
((x>200)*(x<=850)*286)+
((x>850)*(x<=900)*(286+2*(x-850)/50))+
((x>900)*(x<=1000)*(288+4*(x-900)/100))+
((x>1000)*(292+20*((x-1000)/3000)));
qt[]=(0.0025*(x<=900))+
((x>900)*(x<=1000)*(0.0025-(0.002*(x-900)/100)))+
((x>1000)*(x<=2000)*(0.0005+(0.0025*(x-1000)/1000)))+
((x>2000)*(x<=3500)*(0.003-(0.001*(x-2000)/1500)))+
((x>3500)*(0.002-(0.0005*(x-3500)/500)));
}
}
}
}
event Diffusion(i++){
boundary({u,v,T1,qt});
nn=0;
scalar rx[],ry[],rT1[],T1f[],rqt[],b[],w[];
face vector kh[],sqd[],Ri[],fRi[];
double CN,U,bbottom,es,qsl;
kh.x.refine=no_restriction;
sqd.x.refine=no_restriction;
Ri.x.refine=no_restriction;
fRi.x.refine=no_restriction;
foreach()
b[]=((T1[]-Tref)*(9.81/Tref))*(1+(0.608*qt[]));
foreach(){
w[]=-0.005*(((x/1000)*(x<=1000))+(x>1000))*(((t/3600)+16-24)>=24);
rx[]=0.000139*(v[]-Vgeo);
rx[]-=w[]*((u[1]-u[-1])/(2*Delta));
rT1[]=-w[]*((T1[1]-T1[-1])/(2*Delta));
ry[]=0.000139*(Ugeo-u[]);
ry[]-=w[]*((v[1]-v[-1])/(2*Delta));
rqt[]=-w[]*((qt[1]-qt[-1])/(2*Delta));
if (x<Delta){ //Compute Surface fluxes
bbottom=((T1bottom-Tref)*(9.81/Tref))*(1+(0.608*qt[]));
TK1= T1bottom*exner;
es=610.78*exp(17.27*(TK1-273.16)/(TK1-35.86));
qsl= (Rd/Rv)*(es/(97200.-(1-((1-Rd/Rv)*es))));
qtsat=qsl*((1+((sq(L)/(Rv*cpd*sq(TK1)))*qt[]))/(1+((sq(L)/(Rv*cpd*sq(TK1)))*qsl)));
CN=sq(0.4/log((x)/zo));
if (b[]>bbottom){ //Stable
Cm=CN*fris(((x-zo)*(b[]-(bbottom))/(sq(u[])+sq(v[]))));
Ch=Cm;
Cq=Cm*0.025;
}
else{ //Unstable
CN = sq(0.4/log((x)/zo));
Cm = CN*friubm((x-zo)*(b[]-(bbottom))/(sq(u[])+sq(v[])),CN);
Ch = CN*friubh((x-zo)*(b[]-(bbottom))/(sq(u[])+sq(v[])),CN);
Cq = Ch*0.025;
}
U=sqrt(sq(u[])+sq(v[]));
rx[]-=(u[]*Cm*U)/Delta;
ry[]-=(v[]*Cm*U)/Delta;
rT1[]-=((T1[]-T1bottom)*Cm*U)/Delta;
rqt[]-=((qt[]-qtsat)*Cq*U)/Delta;
}
}
boundary(all);
foreach_face()//Compute turbulent diffusivities
{
sqd.x[]=(sq((u[]-u[-1])/(Delta))+sq((v[]-v[-1])/(Delta)));
Ri.x[]= ((b[]-b[-1])/(Delta))/(sqd.x[]+0.00001);
if (Ri.x[]<0)
fRi.x[]=friu(Ri.x[]);
else
fRi.x[]=fris(Ri.x[]);
kh.x[]=sq(min(0.4*x,lz))*(sqrt(sqd.x[]))*fRi.x[];
}
boundary({kh.x});
dt=dtnext(DT);
// Time integrate and log the total number of Multigrid Cycles
mgb=diffusion(u,dt,kh,rx);
nn+=mgb.i;
mgb=diffusion(v,dt,kh,ry);
nn+=mgb.i;
mgb=diffusion(T1,dt,kh,rT1);
nn+=mgb.i;
mgb=diffusion(qt,dt,kh,rqt);
nn+=mgb.i;
}
output
We output the profiles at the specified time.
event outprof (t = 79200){
char fname[99];
sprintf(fname, "%d", j);
FILE * fp = fopen (fname, "w");
foreach()
fprintf(fp, "%g\t%g\n", x,T1[]);
fclose(fp);
}
Stop
The simulation is stopped right after we have outputted the desired profiles.
event stop (t = 79202){
return 1;
}
event adapt(i++;t<He*3600){
if (j!=3)
adapt_wavelet({u,v,T1},(double[]){eu,eu,eb},maxlevel,2,{u,v,T1,qt});
}
Results
The profiles:
set yr [0:1100]
set xr [286:293]
set ylabel 'Height [m]'
set xlabel 'Temp [K]'
set key box bottom right
plot '0' u 2:1 w l lw 2 t 'Default Lz', \
'1' u 2:1 w l lw 2 t 'Decreased Lz', \
'2' u 2:1 w l lw 2 t 'Increased Lz', \
'3' u 2:1 w l lw 2 t 'Default Lz + equidistant grid'
The plot looks like a less-cared for version of the one presented in the rebuttle to ref. 1 and 2.