FEHM Version 3.4.0 - September 2019

The major changes in the FEHM application from V3.3.1 to this release V3.4.0 are: salt and GDKM development, various modifications to improve the performance of FEHM simulation, fixing several known code bugs, and salt_test in V&V test suite was updated with improved theoretical solution (theoretical_result.txt).

See PDF V3.4 Release Notes for additional discussions:

BOUNDARY CONDITIONS FOR THE AIR-WATER-HEAT MODULE IN FEHM AND A DESCRIPTION OF VERIFICATION PROBLEMS
DEVELOPMENT OF A HIGH TEMPERATURE TABULAR WATER PROPERTIES CAPABILITY IN FEHM

Enhancements

Nonlinear rock heat capacity including rock melting. This capability allows the rock heat capacity with temperature dependence to be used in simulations. The numerical formulation for rock energy was modified to allow internal energy changes that represent rock melting. Several different heat capacity/melting models are available including input in tabular form. This capability will be fully tested and input/output described when the High Temperature version of FEHM is merged in the next release.
Negative volume. A simple check for negative volumes associated with inconsistent local element node ordering was implemented.
Humidity boundary conditions. Humidity boundary conditions were added to the BOUN macro.
Gravity term correction. There are instances in multiphase simulations where a single phase gridblock is next to a multiphase block. In this case, it is possible that the density of one of the phases is not calculated for the single phase block. This affects the gravity term and is now corrected. (Seems to save a few iterations.) First discovered in TOUGH code at the University of Auckland.
Salt project. Modified several of the salt routines. These changes were primarily for improvements to output. There are some computational efficiency improvements (saltctr).
SRC_CH2017 (general). Moved vapor pressure lowering calls and improved the free drainage boundary condition for isothermal flow.
Principal stress output. Changes the output for principal stress. Moved some arrays to use modules for improved efficiency.
Saved Zone capability. This allows zones to be saved once and for all. This allows sets of zones not to lose nodes because of overlapping nodes defined sequentially. Useful for setting stratigraphy and boundary conditions and for outputting separate contour zones.
GDKM. This is a major update to the GDKM capability. Corrected known bugs, simplified input, and added models. The models include the generic fracture model (similar to dpdp) and three directional fracture models.

Code Bugs Fixed

varchk_2.f - fixed a bug in the variable update portion of routine for AWH and SC phase
resetv.f - fixed a bug. Now allows a timestep restart when a SC phase is present
Residual.f – bug fix so that 135 Xe concentration test passes
Sub_gmres1.f - bug fix so that 135 Xe concentration test passes
thrmwc.f- continued improvement of ngas physics (AWH) boundary conditions, debugged and modified high temperature table
vcon.f - bug fix for earlier implementation of thermal conductivity for porosity above 0.3873 to FEHM. Included temperature dependence calculation and accurate linear slope equation for porosity above 0.3873 situation.

Known Issues

08/29/2019: problems compiling FEHM with Intel Visual Fortran Checking out the flag “Check Routine Interfaces” in the “Diagnostics” option of the Fortran workspace properties can let compilation of FEHM go through. There are two types of issues: 1) related to the setbit subroutine - the dummy arguments are all integer (4 bytes) while the third argument in the calling is of type real*8. Even if compilation goes through, it is not sure that results could be reliable; 2) there are many cases where the size of dummy arrays are larger than the ones in the calling sequence - but this is not a problem.

Code Changes

------------ Code changes 10/07/2017:
  
GRAV_TERM:
comji.f - added array grav_wgt(:)
flow_boundary_conditions.f - fixed error on specified air pressure
geneqc.f - set grav wgt = 1 if one of the densities (i.e. rovf(i) le.0)
varchk.f - adjustment to default saturation tolerance
  
SALT_OLivella:
comai.f - added global parameter 'iaprf' to manage Olivella intrinsic permeability reduction function; added global para meters in support of code modifications
comdi.f - added globalparameters k0f, bkf, por0ffor the Olivella intrinsic permeability reduction function
data.f - zeroed parameter 'iaprf'
porosi.f - Added coding to enable call to Olivella permeability function (including perm initialization); Corrected the number of loop calls and improved cpu.
saltctr.f - Added input coding for Olivella thermal conductivity (and warning in model 4 or 5 are not used with salt). Also added input for Olivella intrinsic perm function. Appended subroutine perm_olivella to the end of saltctr.f. Added conditional to Olivella perm function to enable calculation of a single gridblock to be consistent with loop on porosi.
fehmn_pcx.f - changed comment only near call to porosi()
flow_boundary_conditions.f - added code so that a specified saturation entered in the boun macro works with Richards' equation
vcon.f - Corrected existing model to implement Olivella thermal conductivity function.
(still some question as to Temp units K or C (Olivella, 2011 seems to be C)
Can't find a reference temperature.
  
PRINCIPAL_STRESS:
write_avs_head_s.f - changed printout order of header for principal stresses
SRC_CH2017:
psatl.f - modified vapor pressure to use Sparrow for salt repository
vaporl.f - added coding to include sparrows 2003 equation
thrair.f - modified free drainage BC
SV_ZONE:
zone.f - added save zone capability. Can define zone and then use them later in the simulation.
They won't disappear. Also supports special SoilVision output.
write_avs_node_v.f - added coding to support SoilVision output
write_avs_node_s.f - added output for gdkm with blanking and no dual grid
write_avs_node_con.f - added output for gdkm with blanking and no dual grid. Added mods for SoilVision support.

SRC_CH2017:
psatl.f - modified vapor pressure to use Sparrow for salt repository
vaporl.f - added coding to include sparrows 2003 equation
thrair.f - modified free drainage BC
  
SV_ZONE:
zone.f - added save zone capability. Can define zone and then use them later in the simulation.
They won't disappear. Also supports special SoilVision output.
write_avs_node_v.f - added coding to support SoilVision output
write_avs_node_s.f - added output for gdkm with blanking and no dual grid
write_avs_node_con.f - added output for gdkm with blanking and no dual grid. Added mods for SoilVision support.

read_avs_io.f - enable SOIL VISION output
initdata2.f - modified to use saved zones.
fehmn_pcx.f - perform restart if equation normalization fails
allocmem.f - added memory allocation for new features and also for SoilVision arrays.
user_ymp.f - moved array elem_temp(:) to use module combi.f

NONLIN_CP:
comdi.f - added variables associated with rock melting
enthp.f - modified for heat conduction only. (might have to set temperature with hflx macro)
input.f - added input call to vrock_ctr.f for rock melting
thermw.f - added coding for rock melting
thrmwc.f - added coding for rock melting. Added coding to implement new humidity related boundary conditions.
vrock_ctr,f - new controller for rock melting
wrtout.f - corrected ethalpy output in *.out file for heat conduction
inrock.f - changed default values to real*8 (i.e. 2500. to 2500.d0)
  
HUMIDITY:
New humidity BC development:
Note: most of the files were changed to support developments in humidity initialization, output and boundary conditions. 'humidity' is used for the relative humidity fraction.
model_setup.f - added boun sub keywords 'huf' (flowing humidity),'hu' (fixed humidity),
th (humidity temperature), ph (humidity pressure)
co2ctr.f - set relative humidity fraction at input temperature for gas phase, print out relative humidity fraction
flow_humidity_bc.f (new) - calculates the air and water mass fractions of an inflowing gas flow so that the inflow has a specified relative humidity.
wrtcon.f - added coding to output more detail on SIA iterations
flow_boun.f - added coding for new humidity boundary conditions.

GDKM_MODS:
add_gdpm.f - added input for newer gdkm directional models
combi.f - added array for gdkm modifications
coneq1.f - modified diffision/dispersion to have correct weighting with gdkm models
check_rlp.f - made temp array sizes consistent with corresponding variable arrays
dvacalc.f - move dva_save(:) to use module comci.f, Added GDKM weighting
comci.f - added array dva_save(:)
gdkm_connect.f - modified volume fraction calcs
gdkm_volume_fraction_interface.f - used volume fractions to calculate gdkm inteface weights
geneq1.f - modified for new gdkm models, modified grav term weighting for zero phase density
geneq2.f - modified for new gdkm models, no gav weighting mods for isothermal
ingdpm.f - modified input for new gdkm models, major mofification for gdkm input
scanin.f - scan input for new gdkm input. Added coding to identify new humidity related keywords in the boun macro.
startup.f - added and rearranged calls to gdkm calcs
avsio.f - added parameter iogdkmblank for gdkm graphics
setparams.f - added more checks for gdpm and gdkm problems
innode.f - added more checks on input when gdkm not used
solve_dual.f - changes some allocate/deallocate statements to clean up subroutine call
csolve.f - corrected do loops for variable gdkm node number, modified to printout SIA iterations

NEG_VOL:
anonp.f - Added check and correction for negative volumes associated with inconsistent local element node ordering.
gncf3.f - Added check and correction for negative volumes associated with inconsistent local element node ordering.
-------- Code changes 10/29/2017:
varchk.f - adjustment to default saturation tolerance. For pure water, calculated T for 2-phase conditions where s is the variable. Needed for rock melting.
thermw.f - added coding for rock melting. Moved where rock internal energy is calculated.
(so T is available for 2-phase problems)
-------- Code changes 12/13/2017:
bnswer.f - small change for managing "out of bounds restart" parameter (mlz)
cntlin.f - added additional unit numbers for hi temp water data table (like co2 tables)
comai.f - added parameters for tracking number of supercritical water nodes
comsi.f - added variable arrays for tabular mechanical properties
comxi.f - increase size of arrays that identify input tabular data (support for hi temp water data)
dated.f - changed the name of the FEHM version to "FEHM V3.3.J HT 17-12-12" (just a suggestion)
enthp.f - changed to compatibility with High Temp code
psatl.f - added derivative wrt pressure for saturation temperature calculation
fehmn_pcx.f - added another use module (use property_interpolate_1) for hi temp water. Added calls to manage "saved zones".
h2o_properties_new.f - New routine to manage hi temp water properties. Based on Pawar's co2_properties.f
inpres.f - added code to allow super critical water phase state
inrlp.f90 - only a few comments but no active line changes
interpolate_2.f90 - New routine for hi temp water property interpolation. Based on John Doherty's routine for CO2 property interpolation.
iofile.f - added coding to manage opening of hi temp water interpolation table
outbnd.f - corrected multiplier for permeability units. (outputted value only)
psatl.f - added coding to include hi temp water table interpolation
rlperm.f - added coding for supercritical water phase (1 line)
scanin.f - initialized iwater_table = 0, determine if using hi temp table
sther,f - set larger T and P bounds if using water interpolation table
stress_mech_props.f - added coding for tabular input of mechanical properties
stressctr.f - added coding for tabular mechanical properties
thermw.f - added coding for hi temp tabular water properties and temperature dependent
(including rock melting) rock/soil properties
thrmwc.f - added coding for hi temp tabular water properties and temperature dependent
(including rock melting) rock/soil properties
varchk.f - checked for supercritical water phase and phase change
writeio.f – added: use comco2, only: icarb
wrtout.f - added output for supercritical water phase information
Notes:
1. The change to rlperm.f was simple: if (ieosd.eq.4) ieosd =1. Not sure, but this may be needed in routines that are related to the keyword 'rlpm'
2. Parameter "iwater_table" refers to hi temp tabular water properties.

--------- Code changes 02/10/2018:
add_gdpm.f- modified geomtriv gpkm terms to be closer to dpdp
geneq1.f - Modified connection between primary and secondary nodes to solve the heat conduction GDKM problem
geneq2.f - general cleanup of primary and secondary connection for consistency
geneqc.f - general cleanup of primary and secondary connection for consistency
write_avs_node_s.f - Fixed overflow of gdkm output for infiltration_gdkm verification problem - now passes the tests
-------- Code changes 03/19/2018:
zone.f - changed definition of zone keyword 'all' from setting every node (including secondary porosity nodes) to the designated zone to just the primary nodes. The secondary node are set '+ 100' or '+ 1000' if primary zone number > 100
infiles.f - modified so the zone file identified in the control file correctly sets the secondary porosity nodes.
ingdpm.f - slight modification to insure old style gdkm input always matches the 'dpdp' non directional fracture type
add_gdpm.f - corrected error to insure the correct area term is used in the non-directional fracture model (gdkm_dir = 0) when 'gdkm new' is used
-------- Code changes 06/20/2018:
geneq2_uz_wt.f - added GDKM coding
write_avs_node_s.f - added blanking for more variables when GDKM used
add_gdpm.f - tweaked some geometric factors for directional fracture models
geneq1.f - changed permeability weighting for primary-secondary weighting for GDKM nodes
geneqc.f - changed permeability weighting for primary-secondary weighting for GDKM nodes
coneq1.f - small GDKM changes
setparams.f - initialized some GDKM parameters
geneq2.f - changed permeability weighting for primary-secondary weighting for GDKM nodes
ingdpm.f - changed some input for new gdkm models
gdkm_volume_fraction_interface.f - corrected some minor errors in GDKM interface weighting.
-------- Code changes 09/09/18:
fehmn_pcx.f - moved the following from startup
if (ianpe.eq.0) then
if (irun.eq.1.and.inobr.eq.0) call sx_combine(1)
  else
if (irun.eq.1.and.inobr.eq.0) call sx_combine_ani(1)
endif
so that beaking connections between sources or sinks is done consist with boun macro and the flow macro
startup.f - moved the same lines above from startup to fehmn_pcx
model_setup.f - added boun keyword 'fen' for flowing enthalpy. Definition is analogous to flowing temperature
scanin.f - added search for boun keyword 'fen'
zone.f - corrected some coding when creating 'saved' zones; not in common usage yet
---------- Code changes 09/25/18:
thrmwc.f - Added coding so that if, on a gridblock with a specified total pressure, the air partial pressure becomes greater than the specified total pressure is released.
Significant change for ngas problems.
varchk.f - phase_mult_hma : changed from 1.15 to 1.15
comai.f - added a few global variables
rlperm.f - added coding so that macro rlp can be called multiple times
scanin.f - added coding to support multiple rlp calls.
initdata2.f - added coding (padding) to support multiple rlp calls
data.f - add and initialize some variables associated with multiple rlp calls
setparams.f - initialized counter for nultiple rlp calls
-------- Code changes 10/26/18:
zone.f - corrected mismatched header lines in calls to zone_elem (a subroutine in zone.f); separated arrays in deallocate statement to check in allocated; changed zonesavename to character*30
input.f - corrected header line in call to zone_elem; changed zonesavename to character*30
initdata2.f - changed zonesavename to character*30
write_avs_node_con.f - changed zonesavename to character*30
write_avs_node_s.f - changed zonesavename to character*30
write_avs_node_v.f - changed zonesavename to character*30
saltctr.f - changed variable in initdata2 call to include macroread(23) in header; added line 'macroread(23) = .false.' before call to initdata2
comai_1.f - only added comment line to identify macroread(23) as associated with perm_olivella
-------- Code changes 11/01/2018 – 11/27/2018:
thrmwc_2.f - added coding to support new uses of boun keywords and improve numerical performance in the "source term code" section
flow_boun.f - added and used array qaxf(i) to store boun air fraction data
flow_boundary_conditions.f - Implemented array qafx() for applying air fraction boundary conditions. (boun keyword xfa)
model_setup.f - modified array usage for 'fxa' boun keyword
comdi.f - added an allocatable array sk_temp used in thrmwc.f; added an allocatable array qaxf used in flow_boun.f.
comai_1.f - added comment for macroread(24) for macro den
saltctr.f - added macroread(23) = .true. to indicate perm_olive read
hydiss.f - changed macroread(7) to macroread(25)
den_vis_spatial.f - changed macroread(12) to macroread(24)
--------- Code changes 11/30/2018:
flow_boundary_conditions.f - fixed logic that allowed access to an unallocated array when fixed temperatures and fixed enthalpy flows occur in the same simulation.
--------- Code changes 12/10/2018:
thrmwc.f- continued improvement of ngas physics (AWH) boundary conditions
debugged and modified high temperature table
varchk.f - improved phase change code for ngas physics with high temperature table
--------- Code changes 12/30/2018:
thrmwc_5.f - continued improvement of ngas physics (AWH) boundary conditions- fixed some more problems with the supercritical water phase
varchk_2.f - fixed a bug in the varible update portion of routine for AWH and SC phase
resetv.f - fixed a bug, now allows a timestep restart when a SC phase is present.
air_cp.f - attached a subroutine air_sol to air_cp to manage solubility models of air in liquid water. For now limited to the original Henry's law.
--------- Code changes 01/12/2019 (merge with Peter Johnson’s development):
allocmem.f - OK to merge, extra variables only
cappr.f - OK to merge, extra code for main PJ stuff fixed problem with ifblock for cap model 1
comdi.f- merged newer code with PJ's now identical in src after 102518
rlperm_1a.f - added PJ code (extra relperm models)
saltctr.f - OK to merge, minor extra PJ code
vcon.f - OK to merge, minor extra PJ code
air_cp.f - added maximum liquid air mass fraction (xnl = 0.1)
varchk_2a.f - small change to phase change criteria
--------- Code changes 01/12/2019:
air_cp.f - corrected henry's law constant
vcon.f - corrected model 4 (PJ merge)
thrmwc_6.f - more tweaks of derivatives for ngas physics
--------- Code changes 01/28/2019:
model_setup.f - changed coding to distinguish between dsa and fxa keywords when multiple boun models exist with sa and fxa respectively
flow_boun.f - added more models that allow the use of ti_linear and corrected an error with distributed source air (dsa). keywords allowed with ti_linear: fd, sw, dsw, dsco2, sa, dsa, se, dse, t
--------- Code changes 02/17/2019:
flow_boun.f - added more models that allow the use of ti_linear and corrected an error with distributed source air (dsa). keywords allowed with ti_linear: fd, sw, dsw, dsco2, sa, dsa, se, dse, t, ft, hd, hdo, s
--------- Code changes 03/01/2019 – 05/22/2019:
startup_1.f-no substantive changes (a debug line)
write_avs_node_s.f - added threshhold for saturation
thrair_1.f - modified seepage face boundary condition (flo3 input) for use with Richards equation solution
steady.f - modified tolerance changing algorithm when using flowrate to determine steady state
comsteady.f - added global variables to support changes made to steady.f
user_ymp1.f - added some operations for manipulating zones. (these do not affect any tests)
fehmn_pcx_1.f - comment lines only
comai_1.f - added some variables
incoord.f - small change to indicate bad grid
wrtout_1.f - added output (additional columns) for zone numbers and phase state for water. Improved accuracy of phase change count.
flxz_2.f - added output (additional column) for average zone saturation in output for water
gradctr_1.f - corrected error in visualization when grad macro is used. Does not affect simulation but only zones in some contour output files when grad macro is used
--------- Code changes 05/22/2019 – 06/11/2019:
user_ymp.f - Changed the general 'a' format, allowed in Intel Fortran, to a format that specifies the number of characters (ie a8) so that the routine compiles in linux.
--------- Code changes 06/13/2019:
residual.f - line 187, anorm_tol change to 1.d-90
sub_gmres1.f - line 370 change tols = 1.d-30, line 427 change 1.d-20 to 1.d-60
--------- Code changes 06/27/2019:
wrtout_1.f - added output (additional columns) for zone numbers and phase state for water. Improved accuracy of phase change count. Added output for zone and phase state information for printout nodes
flxz_2.f - added output (additional column) for average zone saturation in output for water. Corrected memory error when idof = 13 (s always=1)
setconnarray_1.f - added code that writes a file (red_fac_itfc.txt) for visualization of the interface nodes.
--------- Code changes 07/08/2019:
combi.f - added global variables for the number of interface connections affected by itfc macro.
scanin.f - saved the actual number of classes (lines) of interface connections values entered on input for itfc macro.
setconnarray.f - fixed error in colloid part of itfc. The code was previously going through the colloid loop even if there was no colloid input. Added vtk output along with text file for visualizing the interface positions for the reduction factors.
---------- Code changes 08/06/2019:
flxz_2.f - added output (additional column) for average zone saturation in output for water. Corrected memory error when idof = 13 (s always=1)
airctr_2.f - modified phase change criteria 2-phase to gas
combi_1.f - added variables associated with itfc macro (visualization)
gradctr_1.f - correction for case of reading gradient info from aux file.
(As happens when gradctr used multiple times.)
pest_1.f - added for information for flowrate columns
scanin_1.f - added logical variables associated with itfc macro
setconnarray_3.f - added option to find all connections to a particular zone for setting interface multipliers
write_avs_node_con_2.f - significant changes for transport species to get heading correct for tec files. (This could cause some problems.)
wrtout_1.f - write out additional columns with zone number and phase state (minor tweaks over last version)
reset_1.f - reset saturation = 0 ti ieos=2 for Richard’s Eq only (could be significant)
steady_1.f - minor improvements to steady macro. Forced variable tolerance to be monotonic.
comsteady_1.f - added variable which saves last timestep tolerance
comdi_1.f - added last timestep source array
thrair_2.f - added capability to manage maximum source changes for constant pressure nodes and seepage faces (currently disabled).
---------- Code changes 08/28/2019:
vcon.f - bug fix for pjjohnson's implementation of thermal conductivity for porosity above 0.3873 to FEHM. Included temperature dependence calculation and accurate linear slope equation for porosity above 0.3873 situation.
---------- Code changes 09/03/2019:
airctr.f - only comments
comai.f - changed sat_ich from integer to real*8 (sat value below which head is set to zero for output)
comdi.f - added allocatable arrays for renubering algorithms
comsteady.f - added parameters related for improved algorithm for changing eq tolerance
contr.f - corrected errors when using macro 'chea' for output in head for non-head problems
convctr.f - fixed conflict with with crl(1,1) usage for density in isothermal problems and output using subroutine water_density.
diagnostics.f - added zone id to timestep restart information
fehmn_pcx.f - eliminated negative porosity node (eliminated nodes) from phase change count
headctr.f - implimented a saturation id for head output when "chea" macro used
input.f - corrected "chea" default saturation for zero head designation
pest.f - added output for zone id, average sat of a zone, and velocity
renum.f - added infrastructure for some renumbering algorithms
reset.f - if richard's eq enabled, only allow reset to liquid or two-phase condition (gas phase included in two phase)
startup.f - added some comments
steady.f - improved the algorithm for reducing eq tolerance as flow balance converges. Added input parameters:
"tstr" : specify time(days) to start new tolerance reduction algorithm
"sacc" : set tolerance for nodal mass change
thrair.f - added some comments

V3.4.0 Release Notes

FEHM Version 3.4.0 - September 2019

Enhancements

Code Bugs Fixed

Known Issues

Code Changes