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;

    1. Default run
    2. Using a decreased max. mixing length lz
    3. Using an increased max. mixing length lz
    4. 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'
    (script)

    (script)

    The plot looks like a less-cared for version of the one presented in the rebuttle to ref. 1 and 2.