programmer's documentation
Data setting for the groundwater flow module

Introduction

The Hydrogeology module of is a numerical model for water flow and solute transport in continuous porous media. The flow part is based on the Richards equation, derived from the Darcy law and the conservation of mass. The transport part is based on the the classical advection-diffusion equation, slightly modified to account the specificities of underground transport.

This module can be used to simulate transferts of water and solutes in several saturated and/or unsaturated porous media. The flow part can be steady or unsteady, with scalar or tensorial permeabilities and allows any type of soil water retention model, such as the van Genuchten model. The transport part considers dispersion, sorption and radioactive decay.

Physical concepts and equations are presented in the theory guide.

The groundwater flow module is recent, and thus, has few limitations:

Activation of the module

The module can be activated in the usppmo routine in cs_user_parameters.f90. The corresponding keyword is idarcy in the Module for calculation options module:

! Set groundwater flow module active (1: active, 0: inactive)
ippmod(idarcy) = 1

Note that the activation of the module requires to desactivation the turbulent model in usipph routine in cs_user_parameters.f90 file:

if (ixmlpu.eq.0) then
! For groundwater flow module, turbulent model is set to 0 for security reason
iturb = 0
endif

Specific parameters

When the module is activated, its specific input parameters should be set in the user_darcy_ini1 routine of cs_user_parameters.f90 file. An example is given in cs_user_parameters-richards.f90.

Flow part

The permeability can be isotropic (scalar) or anisotropic (tensor) but all soils will be treated in the same way (isotropic or anisotropic):

! Set permeability to isotropic (0) or anisotropic (1) for all soils
darcy_anisotropic_permeability = 0

The primary variable of the groundwater flow module is the hydraulic head H=h+z. In order to switch easily to the pressure head, the keyword darcy_gravity can be used and the value darcy_gravity_x/y/z will defined the direction:

! Set gravity to pass from H to h. Example for H = h + z:
darcy_gravity = 1
darcy_gravity_x = 0.d0
darcy_gravity_y = 0.d0
darcy_gravity_z = 1.d0

The convergence criteron of the Newton scheme can be set over pressure or over velocity. It is recommended to keep the criteron over pressure:

! Set convergence criteron of the Newton scheme over pressure (0) or over velocity (1).
! It is recommended to keep the criteron over pressure.
darcy_convergence_criterion = 0

Transport part

The dispersion can be isotropic (scalar) or anisotropic (tensor) but all solutes will be treated in the same way (isotropic or anisotropic):

! Set dispersion to isotropic (0) or anisotropic (1) for all solutes
darcy_anisotropic_diffusion = 0

The transient transport can be based on a steady or unsteasy darcian velocity field:

! Set if the transport part is based on a steady (0) or unsteady (1) flow field
darcy_unsteady = 0

Numerical parameters

Specific numerical parameters can be set in usipsu routine of cs_user_parameters.f90 file. An example is given in cs_user_parameters-richards.f90:

Flow part

In the case of soils of very different permeabilities, the resolution of the Richards equation requires a weighted gradient computation for tetrahedral meshes. This option is only available for soils with isotropic permeabilities (darcy_anisotropic_permeability = 0) for now. It as to be coupled with the standard least squares gradient recontruction.

! Set gradient computation to weighted (1) for high permeability ratio in tetrahedral meshes.
! Only works with isotropic permeability and the standard least square gradient reconstruction.
if (darcy_anisotropic_permeability.eq.0) then
iwgrec(ipr) = 1
imrgra = 1
endif

It is recommended to choose low criteria for gradient reconstruction in order to obtain a smooth darcian velocity field for the transport part. For instance:

! Set low criteria for gradient reconstruction to obtain a smooth velocity field
nswrsm(ipr) = 5000
epsrsm(ipr) = 1.d-10
nswrgr(ipr) = 5000
epsrgr(ipr) = 1.d-10
epsilo(ipr) = 1.d-13

Transport part

In the case of soils of very different diffusion (dispersion or molecular diffusion), the resolution of the transport equation requires a weighted gradient computation for tetrahedral meshes. This option is only available for soils with isotropic dispersion (darcy_anisotropic_diffusion = 0) for now. It as to be coupled with the standard least squares gradient recontruction.

if (darcy_anisotropic_diffusion.eq.0) then
do ii = 1, nscal
! Set gradient computation to weighted (1) for high permeability ratio in tetrahedral meshes.
! Only works with isotropic diffusion.
iwgrec(isca(ii)) = 1
! Set higher maximum number of iteration for reconstruction to increase the accuracy
nswrsm(isca(ii)) = 10
enddo
endif

Time parameters

The total number of iterations and the reference time step are also set in this routine.

ntmabs = 500
dtref = 1.d-3

However, the time step can be modified in cs_user_extra_operations.f90 (see cs_user_extra_operations-flux.f90) in order to modif the time step with time:

! Example of time modification
do iel = 1, ncel
dt(iel) = (ttcabs**0.95d0) * 5.d-2
if (dt(iel).lt.5.d-2) dt(iel) = 5.d-2
if (dt(iel).gt.1.d3) dt(iel) = 1.d3
enddo

Physical parameters

Physical parameters can be set in usphyv routine of cs_user_physical_properties.f90 file. This section presents two examples that can be found in in cs_user_physical_properties-richards_sat.f90 and cs_user_physical_properties-richards_unsat.f90.

Note that, in the flow part, depending on the variable darcy_anisotropic_permeability, the permeability storage table is permeability (iostropic case) or tensor_permeability (aniotropic case). For the transport part, the isotropic part of the diffusion (molecular diffusion and isotropic dispersion) is always stored in cpro_vscalt. Only anisotropic dispersion (i.e. darcy_anisotropic_diffusion = 1) is stored in the tensor visten.

Case of two saturated soils

This example shows how to set physical parameters for two fully saturated soils (isotropic or anisotropic permeability) and several solutes with isotropic dispersion:

Flow part

! Set parameters for soil 1
call getcel('SOIL1', ncelt, lstcel)
! Loop on cell of the list
do icelt = 1, ncelt
! Cell number
iel = lstcel(icelt)
! Set saturation and capacity (0 if saturated)
saturation(iel) = 0.6d0
capacity(iel) = 0.d0
! Set permeability
if (darcy_anisotropic_permeability.eq.0) then
permeability(iel) = 1.d-1
else
tensor_permeability(1,iel) = 1.d-1
tensor_permeability(2,iel) = 1.d-1
tensor_permeability(3,iel) = 1.d-2
tensor_permeability(4:6,iel) = 0.d0
endif
enddo
! Set parameters for soil 2
call getcel('SOIL2', ncelt, lstcel)
! Loop on cell of the list
do icelt = 1, ncelt
! Cell number
iel = lstcel(icelt)
! Set saturation and capacity (0 if saturated)
saturation(iel) = 0.4d0
capacity(iel) = 0.d0
! Set permeability
if (darcy_anisotropic_permeability.eq.0) then
permeability(iel) = 5.d-1
else
tensor_permeability(1,iel) = 5.d-1
tensor_permeability(2,iel) = 5.d-1
tensor_permeability(3,iel) = 5.d-2
tensor_permeability(4:6,iel) = 0.d0
endif
enddo

Transport part

For every solute, the isotropic permeability and the delay should be set in all soils. For instance:

! Definition of the isotropic diffusion (dipersion and moleculer diffusion)
call getcel ('SOIL1', ncelt, lstcel)
darcy_isotropic_dispersion = 1.d0
molecular_diffusion = 1.d-6
do icelt = 1, ncelt
iel = lstcel(icelt)
if (ifcvsl.ge.0) then
velocity_norm = sqrt(vel(1,iel)**2+vel(2,iel)**2+vel(3,iel)**2)
cpro_vscalt(iel) = darcy_isotropic_dispersion * velocity_norm + &
saturation(iel) * molecular_diffusion
endif
enddo
! Definition of the isotropic diffusion (dipersion and moleculer diffusion)
call getcel ('SOIL2', ncelt, lstcel)
darcy_isotropic_dispersion = 0.2d0
molecular_diffusion = 1.d-8
do icelt = 1, ncelt
iel = lstcel(icelt)
if (ifcvsl.ge.0) then
velocity_norm = sqrt(vel(1,iel)**2+vel(2,iel)**2+vel(3,iel)**2)
cpro_vscalt(iel) = darcy_isotropic_dispersion * velocity_norm + &
saturation(iel) * molecular_diffusion
endif
enddo
! No delay for both soils
do iel = 1, ncel
delay(iel) = 1.d0
enddo

Unsaturated media

This example shows how to set physical parameters for a single variably saturated soil (isotropic or anisotropic permeability) and a single solute with molecular diffusion, anisotropic dispersivity and sorption. The van Genuchten model, coupled with the Mualem condition, is used to determine the relation between the moisture content and the presure head (h).

Flow part

First the permeability and the van Genuchten parameters are set:

! Set intrinsic permeability (only depends on soil)
if (darcy_anisotropic_permeability.eq.0) then
ki = 1.d0
else
ki_xx = 1.d0
ki_yy = 1.d0
ki_zz = 1.d-1
endif
! Set values of the Van Genuchten model parameters
ks_param = 0.3d0
thetar_param = 0.078d0
thetas_param = 0.3d0
n_param = 1.56d0
m_param = 1-1/n_param !(Mualem condition)
l_param = 0.5d0
alpha_param = 0.036d0

As the van Genuchten law is based on the pressure head (h), the gravity term is added if necessary:

! Switch from hydraulic head (H=h+z) to pressure head (h)
darcy_h = cvar_pr(iel)
if (darcy_gravity.eq.1) then
darcy_h = cvar_pr(iel) - xyzcen(1,iel)*darcy_gravity_x - &
xyzcen(2,iel)*darcy_gravity_y - xyzcen(3,iel)*darcy_gravity_z
endif

In order to compute the capacity, the saturation and the permeability, the saturated and the unsaturated parts are treated differently.

In the saturated part (h>=0), the water content is equal to the saturated water content, the permeability is equal to the saturated permeability and the capacity is equal to zero:

! Saturated part (h<=0)
if (darcy_h.ge.0) then
capacity(iel) = 0.d0
saturation(iel) = thetas_param
if (darcy_anisotropic_permeability.eq.0) then
permeability(iel) = ks_param * ki
else
tensor_permeability(1,iel) = ks_param*ki_xx
tensor_permeability(2,iel) = ks_param*ki_yy
tensor_permeability(3,iel) = ks_param*ki_zz
tensor_permeability(4:6,iel) = 0.d0
endif

In the unsaturated part (h<0), the water content, the capacity and the permeability are function of the pressure head. They are determined thanks to the van Genuchten law:

! Unsaturated part (h<0)
else
tmp_1 = abs(alpha_param*darcy_h)**n_param
tmp_2 = 1.d0 / (1.d0 + tmp_1)
se_param = tmp_2**m_param
capacity(iel) = -m_param*n_param*(tmp_1)* &
(thetas_param-thetar_param)*se_param*tmp_2/darcy_h
saturation(iel) = thetar_param + se_param*(thetas_param-thetar_param)
if (darcy_anisotropic_permeability.eq.0) then
permeability(iel) = ks_param*ki*se_param**l_param*(1.d0-(1.d0-tmp_2)**m_param)**2
else
tensor_permeability(1,iel) = ks_param*ki_xx*se_param**l_param*(1.d0-(1.d0-tmp_2)**m_param)**2
tensor_permeability(2,iel) = ks_param*ki_yy*se_param**l_param*(1.d0-(1.d0-tmp_2)**m_param)**2
tensor_permeability(3,iel) = ks_param*ki_zz*se_param**l_param*(1.d0-(1.d0-tmp_2)**m_param)**2
tensor_permeability(4:6,iel) = 0.d0
endif
endif

Transport part

First, the values of the longitudinal and transversal dispersivity as well as the molecular diffusion of the single solute are set as follow:

! Set values of the longitudinal and transversal dirpersivity
darcy_anisotropic_dispersion_l = 2.d0
darcy_anisotropic_dispersion_t = 1.d-1
! Set value of the molecular diffusion
molecular_diffusion = 1.d-3

The molecular diffusion (isotropic term) is stored in cpro_vscalt and computed as:

! Computation of molecular diffusion of the diffusion term
if (ifcvsl.ge.0) then
call field_get_val_s(ifcvsl, cpro_vscalt)
do iel = 1, ncel
cpro_vscalt(iel) = saturation(iel)*molecular_diffusion
enddo
else
cpro_vscalt => null()
endif

The anisotropic dispersion is stored in visten and computed as:

! Computation of the isotropic dispersivity
do iel = 1, ncel
! Computation of the norm of the velocity
velocity_norm = sqrt(vel(1,iel)**2+vel(2,iel)**2+vel(3,iel)**2)
! Tensorial dispersion is stored in visten
tmp_lt = darcy_anisotropic_dispersion_l-darcy_anisotropic_dispersion_t
visten(1,iel) = darcy_anisotropic_dispersion_t*velocity_norm + tmp_lt*vel(1,iel)**2/(velocity_norm+1.d-15)
visten(2,iel) = darcy_anisotropic_dispersion_t*velocity_norm + tmp_lt*vel(2,iel)**2/(velocity_norm+1.d-15)
visten(3,iel) = darcy_anisotropic_dispersion_t*velocity_norm + tmp_lt*vel(3,iel)**2/(velocity_norm+1.d-15)
visten(4,iel) = tmp_lt*vel(2,iel)*vel(1,iel)/(velocity_norm+1.d-15)
visten(5,iel) = tmp_lt*vel(2,iel)*vel(3,iel)/(velocity_norm+1.d-15)
visten(6,iel) = tmp_lt*vel(3,iel)*vel(1,iel)/(velocity_norm+1.d-15)
enddo

Note that the sorption is considered as a delay with the K_d hypothesis. It is then computed as follow:

! Computation of the sorption as a delay term (K_d hypothesis)
rho = 1.5d0
kd = 1.d-1
do iel = 1, ncel
delay(iel) = 1.d0+rho*kd/saturation(iel)
enddo

Radioactive decay

The radioactive decay is treated as a source term in the transport equation. An example can be found in cs_user_source_terms-richards_decay.f90:

! Set the first order decay coefficient
lambda = 1.d-2
! Set radioactive decay for the first solute
if (isca(iscal).eq.1) then
do iel = 1, ncel
crvimp(iel) = -volume(iel)*lambda
enddo
endif

Initialisation

The initialisation of the variables required for the flow part (hydraulic head H) and transport part (concentration c) can be done globally:

! Global initialisation
do iel = 1, ncel
! Initialisation of the hydraulic pressure H with a contant gradient of 1
! among z axis and -0.01 among x axis
cvar_pr(iel) = 1.d0*xyzcen(3,iel) - 1.d-2*xyzcen(1,iel)
! Null concentration by default
cvar_scal_1(iel) = 0.d0
! Velocity initialisation for security reason
cvar_vel(1,iel) = 0.d0
cvar_vel(2,iel) = 0.d0
cvar_vel(3,iel) = 0.d0
enddo

or by selecting a precise soil:

! Initialisation per group of volume
call getcel('SOURCE', ncelt, lstelt)
do icelt = 1, ncelt
iel = lstelt(icelt)
cvar_scal_1(iel) = 1.d0
enddo
endif

Boundary conditions

For groundwater flows of water and solutes, the undefined type face iindef is used to impose Dirichlet, Neumann and mixte boundary conditions on hydraulic head H (here pressure) and solutes. Several examples can be found in cs_user_boundary_conditions-richards.f90.

Dirichlet boundary conditions

Dirichlet boundary conditions can be used to impose a value for the hydraulic head H and the concentration c at a given boundary face:

call getfbr('FACE1', nlelt, lstelt)
do ilelt = 1, nlelt
! Face number
ifac = lstelt(ilelt)
! Undefined type face
itypfb(ifac) = iindef
! Velocity is not solved but deduced from darcy law, then no BC are required.
! However, no flux BC are imposed for safety.
icodcl(ifac,iu) = 3
icodcl(ifac,iv) = 3
icodcl(ifac,iw) = 3
rcodcl(ifac,iu,3) = 0.d0
rcodcl(ifac,iv,3) = 0.d0
rcodcl(ifac,iw,3) = 0.d0
! Dirichlet BC on hydraulic head (H = h + z) to impose a constant value
icodcl(ifac,ipr) = 1
rcodcl(ifac,ipr,1) = 10.d0
! Dirichlet BC on centration C (g/m^3)
do ii = 1, nscal
icodcl(ifac,isca(ii)) = 1
rcodcl(ifac,isca(ii),1) = 1.d0
enddo
enddo

It can also be used to impose a hydraulic head profile at another face:

call getfbr('FACE2', nlelt, lstelt)
do ilelt = 1, nlelt
! Face number
ifac = lstelt(ilelt)
! Undefined type face
itypfb(ifac) = iindef
! Velocity is not solved but deduced from darcy law, then no BC are required.
! However, no flux BC are imposed for safety.
icodcl(ifac,iu) = 3
icodcl(ifac,iv) = 3
icodcl(ifac,iw) = 3
rcodcl(ifac,iu,3) = 0.d0
rcodcl(ifac,iv,3) = 0.d0
rcodcl(ifac,iw,3) = 0.d0
! Dirichlet BC on hydraulic head (H = h + z) to impose a constant gradient over z axis
! here \f$ \grad H \cdot \vect{z} = -0.1 \f$
icodcl(ifac,ipr) = 1
rcodcl(ifac,ipr,1) = -1.d0 * cdgfbo(3,ifac)
enddo

Neumann boundary conditions

Neumann boundary conditions can be used to impose fluxes at boundaries:

call getfbr('FACE3', nlelt, lstelt)
do ilelt = 1, nlelt
! Face number
ifac = lstelt(ilelt)
! Undefined type face
itypfb(ifac) = iindef
! Velocity is not solved but deduced from darcy law, then no BC are required.
! However, no flux BC are imposed for safety.
icodcl(ifac,iu) = 3
icodcl(ifac,iv) = 3
icodcl(ifac,iw) = 3
rcodcl(ifac,iu,3) = 0.d0
rcodcl(ifac,iv,3) = 0.d0
rcodcl(ifac,iw,3) = 0.d0
! Neumann BC on hydraulic head (H = h + z) to impose a gradient among the surface normal
! Here \f$ \grad H \cdot \vect{n} = 0 \f$
icodcl(ifac,ipr) = 3
rcodcl(ifac,ipr,3) = 0.d0
! Neumann BC on centration C for boundary surface with outward or null normal flow $V_out$:
! It allows to impose a gradient among the surface normal for a diffusive flux as the
! convective flux is defined by the $C*V_out$.
! Here \f$ \grad C \cdot \vect{n} = 0 \f$
do ii = 1, nscal
icodcl(ifac,isca(ii)) = 3
rcodcl(ifac,isca(ii),3) = 0.d0
enddo
enddo

Note that, for the transport part, Neumann boundary conditions can only be used for boundary surface with outward or null normal flow. In both cases, the prescribed flux is the diffusive flux.

Mixte boundary conditions

The mixte boundary conditions (Robin) can be used to impose a concentration flux at an entrance (inward normal flow at a boundary face). The following example explains how to determine the two parameters of the mixte boundary in order to impose a total flux:

call getfbr('FACE4', nlelt, lstelt)
! Pointers to the mass flux
call field_get_key_int(ivarfl(ipr), kbmasf, iflmab)
call field_get_val_s(iflmab, bmasfl)
do ilelt = 1, nlelt
! Face number
ifac = lstelt(ilelt)
! Undefined type face
itypfb(ifac) = iindef
! Velocity is not solved but deduced from darcy law, then no BC are required.
! However, no flux BC are imposed for safety.
icodcl(ifac,iu) = 3
icodcl(ifac,iv) = 3
icodcl(ifac,iw) = 3
rcodcl(ifac,iu,3) = 0.d0
rcodcl(ifac,iv,3) = 0.d0
rcodcl(ifac,iw,3) = 0.d0
! Dirichlet BC on hydraulic head (H = h + z) to impose a constant value
icodcl(ifac,ipr) = 1
rcodcl(ifac,ipr,1) = 10.d0
! In order to impose a radioactive activity R_act (Bq) at a surface S_act (m^2)
! during a period of time T_act, a mixed (or Robin) BC is used.
! To do so, two quantities are required:
! * the velocity at an entrance is Vel = - mass_flux / surf
! or Vel = mass_flux / surf at an exit
! * the reference concentration Cref = R_act / (S_act * dt * Vel)
r_act = 0.5d0
s_act = 50.d0
t_act = 10.d0
! The total flux is imposed from 0 to T_act
vel = - bmasfl(ifac)/surfbn(ifac)
cref = r_act / (s_act * dt(ntcabs) * vel)
if (ttcabs.gt.t_act) cref = 0.d0
do ii = 1, nscal
icodcl(ifac,isca(ii)) = 1
rcodcl(ifac,isca(ii),1) = cref
rcodcl(ifac,isca(ii),2) = vel
enddo
enddo

Flux computations

In order to compute fluxes at innner or boundary surfaces, the file cs_user_extra_operations.f90 can be used. An example can be found in cs_user_extra_operations-richards_flux.f90. It shows how to compute the total fluxes at to different surface and how to write the evolution with time:

pas_iter = 10 ! Number of iterations separating two flux calculations.
if (mod(ntcabs,pas_iter).eq.0) then
iprev = 0
inc = 1
iccocg = 1
eps_geom = 1.d-10
call field_get_key_int(ivarfl(iu), kimasf, iflmas)
call field_get_key_int(ivarfl(iu), kbmasf, iflmab)
call field_get_val_s(iflmas, imasfl)
call field_get_val_s(iflmab, bmasfl)
call field_get_key_int(ivarfl(isca(1)), kivisl, ifcvsl)
call field_get_val_s(ifcvsl, cpro_vscalt)
scalar_diffusion(1:ncel) = cpro_vscalt(1:ncel)
call viscfa (imvisf, scalar_diffusion, viscf, viscb)
! Example of tracer flux computation through an internal surface.
! Fluxes are calculated whithout reconstruction, and with a simple definition of the concentration at faces
! (i.e. mean of the two adjacent concentrations).
! Thus, the fluxes are not coherent with the real fluxes calculated in bilsc2 (variable "flux")
! They are approximated and should not be used to verify the exact conservation of mass
call getfac('PLANE1', nlelt, lstelt)
flux_in = 0.d0
do ilelt = 1, nlelt
ifac = lstelt(ilelt)
iel11 = ifacel(1,ifac)
iel22 = ifacel(2,ifac)
! Test locality of the cell (only for mpi computation)
if (iel22.le.ncel) then
tra_face = 5.d-1*(cvar_scal_1(iel11)+cvar_scal_1(iel22))
tra_diff = cvar_scal_1(iel11) - cvar_scal_1(iel22)
else
! Remove duplicated boundary boundary faces for parallel computation
tra_face = 0.d0
tra_diff = 0.d0
endif
flux_tmp = imasfl(ifac)*tra_face + viscf(ifac)*tra_diff
! We impose a norm on the direction of the flux, based on the direction
! of main coordinates, to be sure that fluxes of all faces are all summed
! in the same direction
if ( (surfac(1,ifac).le.(-eps_geom) ) .or. &
( (abs(surfac(1,ifac)).le.eps_geom).and. &
(surfac(2,ifac).le.(-eps_geom)) ) .or. &
( (abs(surfac(1,ifac)).le.eps_geom).and. &
(abs(surfac(2,ifac)).le.eps_geom).and. &
(surfac(3,ifac).le.(-eps_geom)) ) ) then
flux_tmp = -flux_tmp
endif
flux_in = flux_in + flux_tmp
enddo
! Example of boundary flux computation
! Fluxes are calculated whithout reconstruction, and with a simple definition of the concentration at faces
! (i.e. mean of the two adjacent concentrations).
call getfac('TOP_SOIL1', nlelt, lstelt)
flux_out = 0.d0
do ilelt = 1, nlelt
ifac = lstelt(ilelt)
iel11 = ifacel(1,ifac)
iel22 = ifacel(2,ifac)
! Test locality of the cell
if (iel22.le.ncel) then
tra_face = 5.d-1*(cvar_scal_1(iel11)+cvar_scal_1(iel22))
tra_diff = cvar_scal_1(iel11) - cvar_scal_1(iel22)
else
! Remove duplicated boundary boundary faces for parallel computation
tra_face = 0.d0
tra_diff = 0.d0
endif
!
flux_tmp = imasfl(ifac)*tra_face + viscf(ifac)*tra_diff
! We impose a norm on the direction of the flux, based on the direction
! of main coordinates, to be sure that fluxes of all faces are all summed
! in the same direction
if ( (surfac(1,ifac).le.(-eps_geom) ) .or. &
( (abs(surfac(1,ifac)).le.eps_geom).and. &
(surfac(2,ifac).le.(-eps_geom)) ) .or. &
( (abs(surfac(1,ifac)).le.eps_geom).and. &
(abs(surfac(2,ifac)).le.eps_geom).and. &
(surfac(3,ifac).le.(-eps_geom)) ) ) then
flux_tmp = -flux_tmp
endif
flux_out = flux_out + flux_tmp
enddo
! Compute sum for parallel computations
if (irangp.ge.0) then
call parsom(flux_in)
call parsom(flux_out)
endif
! Write fluxes in file
impout = impusr(2)
if (irangp.le.0) then
open(impout,file='flux_tracer1.dat',position="append")
write(impout,97) ttcabs, flux_in, flux_out
97 format(3g17.9)
close(impout)
endif
endif