Numa3DInitialize.c File Reference

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <basic.h>
#include <arrayfunctions.h>
#include <mathfunctions.h>
#include <boundaryconditions.h>
#include <physicalmodels/numa3d.h>
#include <mpivars.h>
#include <hypar.h>

Go to the source code of this file.

Functions
double	Numa3DComputeCFL (void *, void *, double, double)
int	Numa3DFlux (double *, double *, int, void *, double)
int	Numa3DStiffFlux (double *, double *, int, void *, double)
int	Numa3DSource (double *, double *, void *, void *, double)
int	Numa3DRusanovFlux (double *, double *, double *, double *, double *, double *, int, void *, double)
int	Numa3DRusanovLinearFlux (double *, double *, double *, double *, double *, double *, int, void *, double)
void	Numa3DCalculateStandardAtmosphere_1 (void *, double, double *, double *, double *, double *)
void	Numa3DCalculateStandardAtmosphere_2 (void *, double, double *, double *, double *, double *)
int	Numa3DInitialize (void *s, void *m)

Function Documentation

double Numa3DComputeCFL	(	void *	,
		void *	,
		double	,
		double
	)

Definition at line 9 of file Numa3DComputeCFL.c.

{
  HyPar  *solver = (HyPar*)  s;
  Numa3D *param  = (Numa3D*) solver->physics;

  int     *dim    = solver->dim_local;
  int     ghosts  = solver->ghosts;
  int     ndims   = solver->ndims;
  double  *u      = solver->u;
  int     index[ndims];

  double max_cfl = 0;
  int done = 0; _ArraySetValue_(index,ndims,0);
  while (!done) {
    int p; _ArrayIndex1D_(ndims,dim,index,ghosts,p);
    double drho,uvel,vvel,wvel,dT,rho0,T0,P0,EP,c,zcoord;

    _GetCoordinate_(_ZDIR_,index[_ZDIR_],dim,ghosts,solver->x,zcoord);
    param->StandardAtmosphere(param,zcoord,&EP,&P0,&rho0,&T0);
    _Numa3DGetFlowVars_         ((u+_MODEL_NVARS_*p),drho,uvel,vvel,wvel,dT,rho0);
    _Numa3DComputeSpeedofSound_ (param->gamma,param->R,T0,dT,rho0,drho,EP,c);

    double dxinv, dyinv, dzinv;
    _GetCoordinate_(_XDIR_,index[_XDIR_],dim,ghosts,solver->dxinv,dxinv); /* 1/dx */
    _GetCoordinate_(_YDIR_,index[_YDIR_],dim,ghosts,solver->dxinv,dyinv); /* 1/dy */
    _GetCoordinate_(_ZDIR_,index[_ZDIR_],dim,ghosts,solver->dxinv,dzinv); /* 1/dz */

    double local_cfl[3];
    local_cfl[_XDIR_] = (absolute(uvel)+c)*dt*dxinv; /* local cfl for this grid point (x) */
    local_cfl[_YDIR_] = (absolute(vvel)+c)*dt*dyinv; /* local cfl for this grid point (y) */
    local_cfl[_ZDIR_] = (absolute(wvel)+c)*dt*dzinv; /* local cfl for this grid point (z) */
    if (local_cfl[_XDIR_] > max_cfl) max_cfl = local_cfl[_XDIR_];
    if (local_cfl[_YDIR_] > max_cfl) max_cfl = local_cfl[_YDIR_];
    if (local_cfl[_ZDIR_] > max_cfl) max_cfl = local_cfl[_ZDIR_];

    _ArrayIncrementIndex_(ndims,dim,index,done);
  }

  return(max_cfl);
}

int Numa3DFlux	(	double *	,
		double *	,
		int	,
		void *	,
		double
	)

Definition at line 7 of file Numa3DFlux.c.

{
  HyPar  *solver = (HyPar*)   s;
  Numa3D *param  = (Numa3D*) solver->physics;
  int     i;

  int *dim    = solver->dim_local;
  int ghosts  = solver->ghosts;
  int ndims   = solver->ndims;
  int index[ndims], bounds[ndims], offset[ndims];

  /* set bounds for array index to include ghost points */
  _ArrayCopy1D_(dim,bounds,ndims);
  for (i=0; i<ndims; i++) bounds[i] += 2*ghosts;

  /* set offset such that index is compatible with ghost point arrangement */
  _ArraySetValue_(offset,ndims,-ghosts);

  int done = 0; _ArraySetValue_(index,ndims,0);
  while (!done) {
    int p; _ArrayIndex1DWO_(ndims,dim,index,offset,ghosts,p);
    double drho,uvel,vvel,wvel,dT,dP,rho0,T0,P0,EP,zcoord;

    _GetCoordinate_(_ZDIR_,index[_ZDIR_]-ghosts,dim,ghosts,solver->x,zcoord);
    param->StandardAtmosphere(param,zcoord,&EP,&P0,&rho0,&T0);

    _Numa3DGetFlowVars_     ((u+_MODEL_NVARS_*p),drho,uvel,vvel,wvel,dT,rho0);
    _Numa3DComputePressure_ (param,T0,dT,P0,dP);
    _Numa3DSetFlux_         ((f+_MODEL_NVARS_*p),dir,drho,uvel,vvel,wvel,dT,dP,rho0,T0);

    _ArrayIncrementIndex_(ndims,bounds,index,done);
  }
  
  return(0);
}

int Numa3DStiffFlux	(	double *	,
		double *	,
		int	,
		void *	,
		double
	)

Definition at line 43 of file Numa3DFlux.c.

{
  HyPar  *solver = (HyPar*)   s;
  Numa3D *param  = (Numa3D*) solver->physics;
  int     i;

  int *dim    = solver->dim_local;
  int ghosts  = solver->ghosts;
  int ndims   = solver->ndims;
  int index[ndims], bounds[ndims], offset[ndims];

  /* set bounds for array index to include ghost points */
  _ArrayCopy1D_(dim,bounds,ndims);
  for (i=0; i<ndims; i++) bounds[i] += 2*ghosts;

  /* set offset such that index is compatible with ghost point arrangement */
  _ArraySetValue_(offset,ndims,-ghosts);

  int done = 0; _ArraySetValue_(index,ndims,0);
  while (!done) {
    int p; _ArrayIndex1DWO_(ndims,dim,index,offset,ghosts,p);
    double drho,uvel,vvel,wvel,dT,dP,rho0,T0,P0,EP,zcoord;

    _GetCoordinate_(_ZDIR_,index[_ZDIR_]-ghosts,dim,ghosts,solver->x,zcoord);
    param->StandardAtmosphere(param,zcoord,&EP,&P0,&rho0,&T0);

    _Numa3DGetFlowVars_               ((u+_MODEL_NVARS_*p),drho,uvel,vvel,wvel,dT,rho0);
    _Numa3DComputeLinearizedPressure_ (param,T0,dT,P0,dP);
    _Numa3DSetLinearFlux_             ((f+_MODEL_NVARS_*p),dir,drho,uvel,vvel,wvel,dT,dP,rho0,T0);

    _ArrayIncrementIndex_(ndims,bounds,index,done);
  }
  
  return(0);
}

int Numa3DSource	(	double *	,
		double *	,
		void *	,
		void *	,
		double
	)

Definition at line 7 of file Numa3DSource.c.

{
  HyPar  *solver = (HyPar*)   s;
  Numa3D *param  = (Numa3D*) solver->physics;
  int     i;

  int *dim    = solver->dim_local;
  int ghosts  = solver->ghosts;
  int ndims   = solver->ndims;
  int index[ndims], bounds[ndims], offset[ndims];

  /* set bounds for array index to include ghost points */
  _ArrayCopy1D_(dim,bounds,ndims);
  for (i=0; i<ndims; i++) bounds[i] += 2*ghosts;

  /* set offset such that index is compatible with ghost point arrangement */
  _ArraySetValue_(offset,ndims,-ghosts);

  int done = 0; _ArraySetValue_(index,ndims,0);
  while (!done) {
    int p; _ArrayIndex1DWO_(ndims,dim,index,offset,ghosts,p);
    double drho,uvel,vvel,wvel,dT,rho0,P0,EP,T0,zcoord;

    _GetCoordinate_(_ZDIR_,index[_ZDIR_]-ghosts,dim,ghosts,solver->x,zcoord);
    param->StandardAtmosphere(param,zcoord,&EP,&P0,&rho0,&T0);
    _Numa3DGetFlowVars_((u+_MODEL_NVARS_*p),drho,uvel,vvel,wvel,dT,rho0);
    _Numa3DSetSource_  ((S+_MODEL_NVARS_*p),param,uvel,vvel,drho,rho0);

    _ArrayIncrementIndex_(ndims,bounds,index,done);

    /* some useless statements to avoid compiler warnings */
    dT = wvel;
    wvel = dT;
  }
  
  return(0);
}

int Numa3DRusanovFlux	(	double *	,
		double *	,
		double *	,
		double *	,
		double *	,
		double *	,
		int	,
		void *	,
		double
	)

Definition at line 9 of file Numa3DUpwind.c.

{
  HyPar   *solver = (HyPar*)  s;
  Numa3D  *param  = (Numa3D*) solver->physics;
  int      done;

  int *dim   = solver->dim_local;
  int ghosts = solver->ghosts;

  int bounds_outer[_MODEL_NDIMS_], bounds_inter[_MODEL_NDIMS_];
  _ArrayCopy1D3_(dim,bounds_outer,_MODEL_NDIMS_); bounds_outer[dir] =  1;
  _ArrayCopy1D3_(dim,bounds_inter,_MODEL_NDIMS_); bounds_inter[dir] += 1;

  done = 0; int index_outer[3] = {0,0,0}, index_inter[3];
  while (!done) {
    _ArrayCopy1D3_(index_outer,index_inter,_MODEL_NDIMS_);
    for (index_inter[dir] = 0; index_inter[dir] < bounds_inter[dir]; index_inter[dir]++) {
      int p; _ArrayIndex1D3_(_MODEL_NDIMS_,bounds_inter,index_inter,0,p);
      double udiff[_MODEL_NVARS_],drho,vel[3],dT,rho0,P0,T0,EP,c;
      double zcoordL, zcoordR;

      if (dir == _ZDIR_) {
        _GetCoordinate_(_ZDIR_,(index_inter[_ZDIR_]-1),dim,ghosts,solver->x,zcoordL);
        _GetCoordinate_(_ZDIR_,(index_inter[_ZDIR_]  ),dim,ghosts,solver->x,zcoordR);
      } else {
        _GetCoordinate_(_ZDIR_,(index_inter[_ZDIR_]  ),dim,ghosts,solver->x,zcoordL);
        zcoordR = zcoordL;
      }

      /* Rusanov's upwinding scheme */

      udiff[0] = 0.5 * (uR[_MODEL_NVARS_*p+0] - uL[_MODEL_NVARS_*p+0]);
      udiff[1] = 0.5 * (uR[_MODEL_NVARS_*p+1] - uL[_MODEL_NVARS_*p+1]);
      udiff[2] = 0.5 * (uR[_MODEL_NVARS_*p+2] - uL[_MODEL_NVARS_*p+2]);
      udiff[3] = 0.5 * (uR[_MODEL_NVARS_*p+3] - uL[_MODEL_NVARS_*p+3]);
      udiff[4] = 0.5 * (uR[_MODEL_NVARS_*p+4] - uL[_MODEL_NVARS_*p+4]);

      /* left of the interface */
      param->StandardAtmosphere(param,zcoordL,&EP,&P0,&rho0,&T0);
      _Numa3DGetFlowVars_         ((uL+_MODEL_NVARS_*p),drho,vel[0],vel[1],vel[2],dT,rho0);
      _Numa3DComputeSpeedofSound_ (param->gamma,param->R,T0,dT,rho0,drho,EP,c);
      double alphaL = c + absolute(vel[dir]);

      /* right of the interface */
      param->StandardAtmosphere(param,zcoordR,&EP,&P0,&rho0,&T0);
      _Numa3DGetFlowVars_         ((uR+_MODEL_NVARS_*p),drho,vel[0],vel[1],vel[2],dT,rho0);
      _Numa3DComputeSpeedofSound_ (param->gamma,param->R,T0,dT,rho0,drho,EP,c);
      double alphaR = c + absolute(vel[dir]);

      double alpha = max(alphaL,alphaR);
      fI[_MODEL_NVARS_*p+0] = 0.5*(fL[_MODEL_NVARS_*p+0]+fR[_MODEL_NVARS_*p+0])-alpha*udiff[0];
      fI[_MODEL_NVARS_*p+1] = 0.5*(fL[_MODEL_NVARS_*p+1]+fR[_MODEL_NVARS_*p+1])-alpha*udiff[1];
      fI[_MODEL_NVARS_*p+2] = 0.5*(fL[_MODEL_NVARS_*p+2]+fR[_MODEL_NVARS_*p+2])-alpha*udiff[2];
      fI[_MODEL_NVARS_*p+3] = 0.5*(fL[_MODEL_NVARS_*p+3]+fR[_MODEL_NVARS_*p+3])-alpha*udiff[3];
      fI[_MODEL_NVARS_*p+4] = 0.5*(fL[_MODEL_NVARS_*p+4]+fR[_MODEL_NVARS_*p+4])-alpha*udiff[4];
    }
    _ArrayIncrementIndex_(_MODEL_NDIMS_,bounds_outer,index_outer,done);
  }

  return(0);
}

int Numa3DRusanovLinearFlux	(	double *	,
		double *	,
		double *	,
		double *	,
		double *	,
		double *	,
		int	,
		void *	,
		double
	)

Definition at line 71 of file Numa3DUpwind.c.

{
  HyPar   *solver = (HyPar*)  s;
  Numa3D  *param  = (Numa3D*) solver->physics;
  int      done;

  int *dim   = solver->dim_local;
  int ghosts = solver->ghosts;

  int bounds_outer[_MODEL_NDIMS_], bounds_inter[_MODEL_NDIMS_];
  _ArrayCopy1D3_(dim,bounds_outer,_MODEL_NDIMS_); bounds_outer[dir] =  1;
  _ArrayCopy1D3_(dim,bounds_inter,_MODEL_NDIMS_); bounds_inter[dir] += 1;

  done = 0; int index_outer[3] = {0,0,0}, index_inter[3];
  while (!done) {
    _ArrayCopy1D3_(index_outer,index_inter,_MODEL_NDIMS_);
    for (index_inter[dir] = 0; index_inter[dir] < bounds_inter[dir]; index_inter[dir]++) {
      int p; _ArrayIndex1D3_(_MODEL_NDIMS_,bounds_inter,index_inter,0,p);
      double udiff[_MODEL_NVARS_],drho,vel[3],dT,rho0,P0,T0,EP,c;
      double zcoordL, zcoordR;

      if (dir == _ZDIR_) {
        _GetCoordinate_(_ZDIR_,(index_inter[_ZDIR_]-1),dim,ghosts,solver->x,zcoordL);
        _GetCoordinate_(_ZDIR_,(index_inter[_ZDIR_]  ),dim,ghosts,solver->x,zcoordR);
      } else {
        _GetCoordinate_(_ZDIR_,(index_inter[_ZDIR_]  ),dim,ghosts,solver->x,zcoordL);
        zcoordR = zcoordL;
      }

      /* Rusanov's upwinding scheme */

      udiff[0] = 0.5 * (uR[_MODEL_NVARS_*p+0] - uL[_MODEL_NVARS_*p+0]);
      udiff[1] = 0.5 * (uR[_MODEL_NVARS_*p+1] - uL[_MODEL_NVARS_*p+1]);
      udiff[2] = 0.5 * (uR[_MODEL_NVARS_*p+2] - uL[_MODEL_NVARS_*p+2]);
      udiff[3] = 0.5 * (uR[_MODEL_NVARS_*p+3] - uL[_MODEL_NVARS_*p+3]);
      udiff[4] = 0.5 * (uR[_MODEL_NVARS_*p+4] - uL[_MODEL_NVARS_*p+4]);

      /* left of the interface */
      param->StandardAtmosphere   (param,zcoordL,&EP,&P0,&rho0,&T0);
      _Numa3DGetFlowVars_         ((uL+_MODEL_NVARS_*p),drho,vel[0],vel[1],vel[2],dT,rho0);
      _Numa3DComputeLinearizedSpeedofSound_ (param->gamma,param->R,T0,rho0,EP,c);
      double alphaL = c + absolute(vel[dir]);

      /* right of the interface */
      param->StandardAtmosphere(param,zcoordR,&EP,&P0,&rho0,&T0);
      _Numa3DGetFlowVars_         ((uR+_MODEL_NVARS_*p),drho,vel[0],vel[1],vel[2],dT,rho0);
      _Numa3DComputeLinearizedSpeedofSound_ (param->gamma,param->R,T0,rho0,EP,c);
      double alphaR = c + absolute(vel[dir]);

      double alpha = max(alphaL,alphaR);
      fI[_MODEL_NVARS_*p+0] = 0.5*(fL[_MODEL_NVARS_*p+0]+fR[_MODEL_NVARS_*p+0])-alpha*udiff[0];
      fI[_MODEL_NVARS_*p+1] = 0.5*(fL[_MODEL_NVARS_*p+1]+fR[_MODEL_NVARS_*p+1])-alpha*udiff[1];
      fI[_MODEL_NVARS_*p+2] = 0.5*(fL[_MODEL_NVARS_*p+2]+fR[_MODEL_NVARS_*p+2])-alpha*udiff[2];
      fI[_MODEL_NVARS_*p+3] = 0.5*(fL[_MODEL_NVARS_*p+3]+fR[_MODEL_NVARS_*p+3])-alpha*udiff[3];
      fI[_MODEL_NVARS_*p+4] = 0.5*(fL[_MODEL_NVARS_*p+4]+fR[_MODEL_NVARS_*p+4])-alpha*udiff[4];

      vel[0] = dT; /* useless statement to avoid compiler warning */
    }
    _ArrayIncrementIndex_(_MODEL_NDIMS_,bounds_outer,index_outer,done);
  }

  return(0);
}

void Numa3DCalculateStandardAtmosphere_1	(	void *	,
		double	,
		double *	,
		double *	,
		double *	,
		double *
	)

Definition at line 5 of file Numa3DFunctions.c.

{
  Numa3D *physics = (Numa3D*) p;

  double R      = physics->R;
  double gamma  = physics->gamma;
  double g      = physics->g;

  /* reference quantities at zero altitude */
  double P0, T0; 
  P0   = physics->Pref;
  T0   = physics->Tref;

  double inv_gamma_m1 = 1.0/(gamma-1.0);
  double Cp = gamma * inv_gamma_m1 * R;

  double theta  = T0;
  *ExnerP = 1.0 - (g/(Cp*theta))*z;
  *P      = P0 * raiseto((*ExnerP),gamma*inv_gamma_m1);
  *rho    = (P0/(R*theta)) * raiseto((*ExnerP),inv_gamma_m1);
  *T      = (*rho) * theta;

  return(0);
}

void Numa3DCalculateStandardAtmosphere_2	(	void *	,
		double	,
		double *	,
		double *	,
		double *	,
		double *
	)

Definition at line 30 of file Numa3DFunctions.c.

{
  Numa3D *physics = (Numa3D*) p;

  double R      = physics->R;
  double gamma  = physics->gamma;
  double g      = physics->g;

  /* reference quantities at zero altitude */
  double P0, T0; 
  P0 = physics->Pref;
  T0 = physics->Tref;

  double BV = 0.01; /* Brunt-Vaisala frequency */
  double inv_gamma_m1 = 1.0/(gamma-1.0);
  double Cp = gamma * inv_gamma_m1 * R;

  double term   = BV*BV*z/g;
  double theta  = T0 * exp(term);
  *ExnerP = 1.0 + (g*g/(Cp*T0*BV*BV)) * (exp(-term) - 1.0);
  *P      = P0 * raiseto((*ExnerP),gamma*inv_gamma_m1);
  *rho    = (P0/(R*theta)) * raiseto((*ExnerP),inv_gamma_m1);
  *T      = (*rho) * theta;

  return(0);
}

int Numa3DInitialize	(	void *	s,
		void *	m
	)

Definition at line 24 of file Numa3DInitialize.c.

{
  HyPar           *solver  = (HyPar*)         s;
  MPIVariables    *mpi     = (MPIVariables*)  m; 
  Numa3D          *physics = (Numa3D*)        solver->physics;
  int             ferr     = 0;

  static int count = 0;

  if (solver->nvars != _MODEL_NVARS_) {
    fprintf(stderr,"Error in Numa3DInitialize(): nvars has to be %d.\n",_MODEL_NVARS_);
    return(1);
  }
  if (solver->ndims != _MODEL_NDIMS_) {
    fprintf(stderr,"Error in Numa3DInitialize(): ndims has to be %d.\n",_MODEL_NDIMS_);
    return(1);
  }

  /* default values */
  physics->gamma  = 1.4; 
  physics->R      = 287.058;        /* J kg^{-1} K^{-1} */
  physics->Omega  = 7.2921150E-05;  /* rad s^{-1}       */
  physics->g      = 9.8;            /* m s^{-2}         */

  physics->Pref   = 101327.0;       /* N m^{-2}         */
  physics->Tref   = 288.15;         /* Kelvin           */

  strcpy(physics->upwind,_RUSANOV_UPWINDING_);

  /* default choice of initial atmosphere */
  physics->init_atmos = 1;

  /* reading physical model specific inputs - all processes */
  if (!mpi->rank) {
    FILE *in;
    if (!count) printf("Reading physical model inputs from file \"physics.inp\".\n");
    in = fopen("physics.inp","r");
    if (!in) printf("Warning: File \"physics.inp\" not found. Using default values.\n");
    else {
      char word[_MAX_STRING_SIZE_];
      ferr = fscanf(in,"%s",word); if (ferr != 1) return(1);
      if (!strcmp(word, "begin")){
          while (strcmp(word, "end")){
              ferr = fscanf(in,"%s",word); if (ferr != 1) return(1);
          if (!strcmp(word, "gamma")) { 
            ferr = fscanf(in,"%lf",&physics->gamma); if (ferr != 1) return(1);
          } else if (!strcmp(word,"R")) {
            ferr = fscanf(in,"%lf",&physics->R); if (ferr != 1) return(1);
          } else if (!strcmp(word,"g")) {
            ferr = fscanf(in,"%lf",&physics->g); if (ferr != 1) return(1);
          } else if (!strcmp(word,"Omega")) {
            ferr = fscanf(in,"%lf",&physics->Omega); if (ferr != 1) return(1);
          } else if (!strcmp(word,"Pref")) {
            ferr = fscanf(in,"%lf",&physics->Pref); if (ferr != 1) return(1);
          } else if (!strcmp(word,"Tref")) {
            ferr = fscanf(in,"%lf",&physics->Tref); if (ferr != 1) return(1);
          } else if (!strcmp(word,"init_atmos")) {
            ferr = fscanf(in,"%d",&physics->init_atmos); if (ferr != 1) return(1);
          } else if (!strcmp(word,"upwinding")) {
            ferr = fscanf(in,"%s",physics->upwind); if (ferr != 1) return(1);
          } else if (strcmp(word,"end")) {
            char useless[_MAX_STRING_SIZE_];
            ferr = fscanf(in,"%s",useless); if (ferr != 1) return(ferr);
            printf("Warning: keyword %s in file \"physics.inp\" with value %s not ",word,useless);
            printf("recognized or extraneous. Ignoring.\n");
          }
        }
        } else {
          fprintf(stderr,"Error: Illegal format in file \"physics.inp\".\n");
        return(1);
        }
    }
    fclose(in);
  }

#ifndef serial
  IERR MPIBroadcast_double    (&physics->gamma      ,1                ,0,&mpi->world); CHECKERR(ierr);
  IERR MPIBroadcast_double    (&physics->R          ,1                ,0,&mpi->world); CHECKERR(ierr);
  IERR MPIBroadcast_double    (&physics->g          ,1                ,0,&mpi->world); CHECKERR(ierr);
  IERR MPIBroadcast_double    (&physics->Omega      ,1                ,0,&mpi->world); CHECKERR(ierr);
  IERR MPIBroadcast_double    (&physics->Pref       ,1                ,0,&mpi->world); CHECKERR(ierr);
  IERR MPIBroadcast_double    (&physics->Tref       ,1                ,0,&mpi->world); CHECKERR(ierr);
  IERR MPIBroadcast_integer   (&physics->init_atmos ,1                ,0,&mpi->world); CHECKERR(ierr);
  IERR MPIBroadcast_character ( physics->upwind     ,_MAX_STRING_SIZE_,0,&mpi->world); CHECKERR(ierr);
#endif

  /* calculate the mean hydrostatic atmosphere as a function of altitude */
  if (physics->init_atmos == 1) {
    physics->StandardAtmosphere = Numa3DCalculateStandardAtmosphere_1;
  } else if (physics->init_atmos == 2) {
    physics->StandardAtmosphere = Numa3DCalculateStandardAtmosphere_2;
  } else {
    if (!mpi->rank) {
      fprintf(stderr,"Error in Numa3DInitialize(): invalid choice of initial atmosphere (init_atmos).\n");
      return(1);
    }
  }
  CHECKERR(ierr);

  /* initializing physical model-specific functions */
  if (!strcmp(solver->SplitHyperbolicFlux,"yes")) 
    solver->dFFunction    = Numa3DStiffFlux;
  else solver->dFFunction = NULL;
  solver->FFunction       = Numa3DFlux;
  solver->ComputeCFL      = Numa3DComputeCFL;
  solver->SFunction       = Numa3DSource;
  if (!strcmp(physics->upwind,_RUSANOV_UPWINDING_)) {
    solver->Upwind        = Numa3DRusanovFlux;
    if (!strcmp(solver->SplitHyperbolicFlux,"yes")) 
      solver->UpwinddF    = Numa3DRusanovLinearFlux;
    else solver->UpwinddF = NULL;
  } else {
    if (!mpi->rank) fprintf(stderr,"Error in Numa3DInitialize(): Invalid choice of upwinding scheme.\n");
    return(1);
  }

  /* set the value of gamma in all the boundary objects */
  int n;
  DomainBoundary  *boundary = (DomainBoundary*) solver->boundary;
  for (n = 0; n < solver->nBoundaryZones; n++)  boundary[n].gamma = physics->gamma;

  count++;
  return(0);
}