//Adapt initial grids

void adaptLESinit(int MAXLEVEL)
{
  unrefine(y>L0/4 && level>=(MAXLEVEL-1),all);
  unrefine(y>200 && level>=(MAXLEVEL-2),all); 
}

//Adapt grid according to solution while running: 
//SGS-flux is compared to resolved-flux to check if the solution is really a LES. If not -> REFINE
//If the SGS-flux is neglegible to resolved flux we are not taking full benefit from LES -> Coarsen
//When fluxes vanish (laminarization) check if Ed.vis < Mol.vis 
//Alse check if there is any resolved flux otherwise i am not sure what to do  


void adaptLESrun(int MAXLEVEL, vector u,scalar T, scalar Evis, int lev[], double yc[]) 
{
  fprintf(ferr,"hoi\n");
  double cur = 0.1, cr=0.5;
  scalar fmr[] , fms[],uh[];	
  foreach()
    {	
      uh[]=sqrt(sq(u.x[])+sq(u.z[]));
      fmr[] = sqrt(sq(u.y[]*(uh[0,1,0]-uh[0,-1,0])/Delta));
      fms[] = sqrt(sq(Evis[]*(uh[0,1,0]-(2*uh[])+uh[0,-1,0])/(sq(Delta))));  	
    }
	
  fprintf(ferr,"hoi\n");
  int k = 0;
  double xp,zp,Url[1<<MAXLEVEL]; 
  int uri = 0;
  while (yc[k] != 0)	
    {
#ifdef _OPENMP
#pragma omp parallel for reduction(+:F) default(shared)
#endif
      double afmr=0 , afms=0 ;
      double yp = yc[k];
      for (int i = 0; i < pow(2,lev[k]); i++)
	{	
	  xp=X0+L0/((pow(2,lev[k]+1))+i*(L0/(pow(2,lev[k]))));	
	  for (int j = 0; j < pow(2,lev[k]); j++)
	    {
	      zp=Z0+L0/((pow(2,lev[k]+1))+i*(L0/(pow(2,lev[k]))));						Point point = locate (xp, yp,zp);
	      afmr += fmr[];
	      afms += fms[];
	    }
	}
      if (cur * afmr > afms )
	{
	  Url[uri] = yc[k];
	  if (uri>0)
	    {	
	      if ( (Url[uri]-Url[uri-1]) == (L0/(pow(2,lev[k]))))
		unrefine(y == ( (Url[uri]+Url[uri]) /2 ),all);
	      
	    }
	  uri++;
	}
      else if (cr * afmr < afms && (afmr/pow(2,2*lev[k])) > 0.01)
	{
	  refine(y>(lev[k]-0.1) && y<(lev[k]+0.1),all);
	}
      k++;
      boundary(all); 
    }
  
}