fix bug on BW eos for density by 1, bypass brine contribution if conc…

…entration field is not specified and 2, give default values to T0, P0 and C0 if not initiated in diamond

CMakeFiles/CMakeOutput.log

@@ -0,0 +1,2 @@
+The system is: Linux - 5.15.0-122-generic - x86_64
+The system is: Linux - 5.15.0-122-generic - x86_64

CMakeFiles/CMakeSystem.cmake

@@ -0,0 +1,15 @@
+set(CMAKE_HOST_SYSTEM "Linux-5.15.0-122-generic")
+set(CMAKE_HOST_SYSTEM_NAME "Linux")
+set(CMAKE_HOST_SYSTEM_VERSION "5.15.0-122-generic")
+set(CMAKE_HOST_SYSTEM_PROCESSOR "x86_64")
+
+
+
+set(CMAKE_SYSTEM "Linux-5.15.0-122-generic")
+set(CMAKE_SYSTEM_NAME "Linux")
+set(CMAKE_SYSTEM_VERSION "5.15.0-122-generic")
+set(CMAKE_SYSTEM_PROCESSOR "x86_64")
+
+set(CMAKE_CROSSCOMPILING "FALSE")
+
+set(CMAKE_SYSTEM_LOADED 1)

ICFERST/Makefile

@@ -0,0 +1,95 @@
+#    Copyright (C) 2006 Imperial College London and others.
+#
+#    Please see the AUTHORS file in the main source directory for a full list
+#    of copyright holders.
+#
+#    Prof. C Pain
+#    Applied Modelling and Computation Group
+#    Department of Earth Science and Engineering
+#    Imperial College London
+#
+#    [email protected]
+#
+#    This library is free software; you can redistribute it and/or
+#    modify it under the terms of the GNU Lesser General Public
+#    License as published by the Free Software Foundation,
+#    version 2.1 of the License.
+#
+#    This library is distributed in the hope that it will be useful,
+#    but WITHOUT ANY WARRANTY; without even the implied warranty of
+#    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+#    Lesser General Public License for more details.
+#
+#    You should have received a copy of the GNU Lesser General Public
+#    License along with this library; if not, write to the Free Software
+#    Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307
+#    USA
+SHELL = /bin/bash
+PACKAGE_NAME = fluidity
+
+FLUIDITY = $(PACKAGE_NAME)
+FC      = mpif90
+FCFLAGS = -ffast-math -frecord-marker=4 -I/usr/include -DHAVE_NUMPY -I/home/rl420/ICFERST1/multifluids_icferst/venv/lib/python3.8/site-packages/numpy/core/include -fno-realloc-lhs -ffree-line-length-none -ffixed-line-length-none  -O3 -I/usr/lib/petscdir/petsc3.12/x86_64-linux-gnu-real/include -I/usr/include/hypre -I/usr/include/suitesparse -I/usr/include/superlu -I/usr/include/superlu-dist -I/usr/include/scotch -I/usr/include/hdf5/openmpi -I/home/rl420/ICFERST2_new_master/multifluids_icferst/include/ -fdefault-real-8 -fdefault-double-8
+
+MPIF90  = mpif90
+MPICC   = mpicc
+MPICXX  = mpicxx
+
+CC	= mpicc
+CFLAGS	= -I/usr/include -I/usr/include/python3.8 -DHAVE_NUMPY -I/home/rl420/ICFERST1/multifluids_icferst/venv/lib/python3.8/site-packages/numpy/core/include -I/usr/lib/petscdir/petsc3.12/x86_64-linux-gnu-real/include -I/usr/include/hypre -I/usr/include/suitesparse -I/usr/include/superlu -I/usr/include/superlu-dist -I/usr/include/scotch -I/usr/include/hdf5/openmpi -I/home/rl420/ICFERST2_new_master/multifluids_icferst/include/ -DHAVE_PETSC  -I/usr/include/vtk-7.1 -I/usr/include/freetype2 -I/usr/include -I/usr/include/x86_64-linux-gnu -I/usr/lib/x86_64-linux-gnu/openmpi/include/openmpi -I/usr/lib/x86_64-linux-gnu/openmpi/include -I/usr/include/python3.8 -I/usr/include/hdf5/openmpi -I/usr/include/jsoncpp -I/usr/include/libxml2 -I/usr/include/tcl   -DHAVE_VTK  -I/usr/include/python3.8 -DHAVE_NUMPY -I/home/rl420/ICFERST1/multifluids_icferst/venv/lib/python3.8/site-packages/numpy/core/include  -O3
+
+CXX	= mpicxx
+CXXFLAGS= -I/usr/include -I/usr/include/python3.8 -DHAVE_NUMPY -I/home/rl420/ICFERST1/multifluids_icferst/venv/lib/python3.8/site-packages/numpy/core/include -I/usr/lib/petscdir/petsc3.12/x86_64-linux-gnu-real/include -I/usr/include/hypre -I/usr/include/suitesparse -I/usr/include/superlu -I/usr/include/superlu-dist -I/usr/include/scotch -I/usr/include/hdf5/openmpi -I/home/rl420/ICFERST2_new_master/multifluids_icferst/include/ -DHAVE_PETSC  -I/usr/include/vtk-7.1 -I/usr/include/freetype2 -I/usr/include -I/usr/include/x86_64-linux-gnu -I/usr/lib/x86_64-linux-gnu/openmpi/include/openmpi -I/usr/lib/x86_64-linux-gnu/openmpi/include -I/usr/include/python3.8 -I/usr/include/hdf5/openmpi -I/usr/include/jsoncpp -I/usr/include/libxml2 -I/usr/include/tcl   -DHAVE_VTK  -I/usr/include/python3.8 -DHAVE_NUMPY -I/home/rl420/ICFERST1/multifluids_icferst/venv/lib/python3.8/site-packages/numpy/core/include  -O3 -I./include
+
+LINKER  = mpicxx  -Xlinker -export-dynamic -Wl,-O1 -Wl,-Bsymbolic-functions   $(CXXFLAGS)
+
+FLLINKER = mpif90  -Xlinker -export-dynamic -Wl,-O1 -Wl,-Bsymbolic-functions   $(FCFLAGS)
+
+MAKE    = make
+AR  = ar
+ARFLAGS = cr
+
+OBJS = $(shell find . -mindepth 2 -maxdepth 2 -name '*.c' -o -name '*.cpp' -o -name '*.[fF]' -o -name '*.F90' | grep -v ^./lib | xargs)
+
+.SUFFIXES: .f90 .F90 .c .cpp .o .a
+
+.f90.o:
+	@echo "        FC $<"
+	$(FC) $(FCFLAGS) -c $<
+.F90.o:
+	@echo "        FC $<"
+	$(FC) $(FCFLAGS) -c $<
+.c.o:
+	@echo "        CC $<"
+	$(CC) $(CFLAGS) -c $<
+.cpp.o:
+	@echo "        CXX $<"
+	$(CXX) $(CXXFLAGS) -c $<
+
+default: lib/libmphase.a  main.o
+
+# cancel implicit rule that tries to make .o from .mod with m2c:
+%.o: %.mod
+
+lib/libmphase.a: $(OBJS)
+	@echo "BUILD libmphase"
+	@echo "    MKDIR lib"
+	@mkdir -p lib
+	@echo "    MAKE src"
+	@cd src; 	 $(MAKE)
+	@$(AR) $(ARFLAGS) lib/libmphase.a src/*.o
+
+clean:
+	@rm -f *.o *.d include/*.mod $(PROG_EXEC) lib/*.a
+	@rm -rf bin/
+	@cd src;   make clean
+	@cd tests; PYTHONPATH=../python ../tools/testharness_ICFERST.py --clean >/dev/null
+	@cd unittest; make clean
+
+
+build_unittest: lib/libmphase.a
+	@cd unittest; $(MAKE)
+
+unittest: build_unittest
+	@echo "RUN bin/tests"
+	@../bin/unittestharness bin/tests

ICFERST/src/Makefile

@@ -0,0 +1,85 @@
+#    Copyright (C) 2006 Imperial College London and others.
+#
+#    Please see the AUTHORS file in the main source directory for a full list
+#    of copyright holders.
+#
+#    Prof. C Pain
+#    Applied Modelling and Computation Group
+#    Department of Earth Science and Engineering
+#    Imperial College London
+#
+#    [email protected]
+#
+#    This library is free software; you can redistribute it and/or
+#    modify it under the terms of the GNU Lesser General Public
+#    License as published by the Free Software Foundation,
+#    version 2.1 of the License.
+#
+#    This library is distributed in the hope that it will be useful,
+#    but WITHOUT ANY WARRANTY; without even the implied warranty of
+#    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+#    Lesser General Public License for more details.
+#
+#    You should have received a copy of the GNU Lesser General Public
+#    License along with this library; if not, write to the Free Software
+#    Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307
+#    USA
+SHELL = /bin/bash
+PACKAGE_NAME = fluidity
+
+FLUIDITY = $(PACKAGE_NAME)
+FC      = mpif90
+FCFLAGS = -I../../include -I../include -J../include -I/usr/include -I/usr/include/python3.8 -DHAVE_NUMPY -I/home/rl420/ICFERST1/multifluids_icferst/venv/lib/python3.8/site-packages/numpy/core/include -I/usr/lib/petscdir/petsc3.12/x86_64-linux-gnu-real/include -I/usr/include/hypre -I/usr/include/suitesparse -I/usr/include/superlu -I/usr/include/superlu-dist -I/usr/include/scotch -I/usr/include/hdf5/openmpi -I/home/rl420/ICFERST2_new_master/multifluids_icferst/include/ -DHAVE_PETSC  -I/usr/include/vtk-7.1 -I/usr/include/freetype2 -I/usr/include -I/usr/include/x86_64-linux-gnu -I/usr/lib/x86_64-linux-gnu/openmpi/include/openmpi -I/usr/lib/x86_64-linux-gnu/openmpi/include -I/usr/include/python3.8 -I/usr/include/hdf5/openmpi -I/usr/include/jsoncpp -I/usr/include/libxml2 -I/usr/include/tcl   -DHAVE_VTK -ffast-math -frecord-marker=4 -I/usr/include -DHAVE_NUMPY -I/home/rl420/ICFERST1/multifluids_icferst/venv/lib/python3.8/site-packages/numpy/core/include -fno-realloc-lhs -ffree-line-length-none -ffixed-line-length-none  -O3 -I/usr/lib/petscdir/petsc3.12/x86_64-linux-gnu-real/include -I/usr/include/hypre -I/usr/include/suitesparse -I/usr/include/superlu -I/usr/include/superlu-dist -I/usr/include/scotch -I/usr/include/hdf5/openmpi -I/home/rl420/ICFERST2_new_master/multifluids_icferst/include/ -fdefault-real-8 -fdefault-double-8 -I/usr/local/include -I./
+
+MPIF90  = mpif90
+
+CC      = mpicc
+CFLAGS  = -I/usr/include -I/usr/include/python3.8 -DHAVE_NUMPY -I/home/rl420/ICFERST1/multifluids_icferst/venv/lib/python3.8/site-packages/numpy/core/include -I/usr/lib/petscdir/petsc3.12/x86_64-linux-gnu-real/include -I/usr/include/hypre -I/usr/include/suitesparse -I/usr/include/superlu -I/usr/include/superlu-dist -I/usr/include/scotch -I/usr/include/hdf5/openmpi -I/home/rl420/ICFERST2_new_master/multifluids_icferst/include/ -DHAVE_PETSC  -I/usr/include/vtk-7.1 -I/usr/include/freetype2 -I/usr/include -I/usr/include/x86_64-linux-gnu -I/usr/lib/x86_64-linux-gnu/openmpi/include/openmpi -I/usr/lib/x86_64-linux-gnu/openmpi/include -I/usr/include/python3.8 -I/usr/include/hdf5/openmpi -I/usr/include/jsoncpp -I/usr/include/libxml2 -I/usr/include/tcl   -DHAVE_VTK  -I/usr/include/python3.8 -DHAVE_NUMPY -I/home/rl420/ICFERST1/multifluids_icferst/venv/lib/python3.8/site-packages/numpy/core/include  -O3 -I../include -I/usr/local/include
+
+CXX     = mpicxx
+CXXFLAGS= -I/usr/include -I/usr/include/python3.8 -DHAVE_NUMPY -I/home/rl420/ICFERST1/multifluids_icferst/venv/lib/python3.8/site-packages/numpy/core/include -I/usr/lib/petscdir/petsc3.12/x86_64-linux-gnu-real/include -I/usr/include/hypre -I/usr/include/suitesparse -I/usr/include/superlu -I/usr/include/superlu-dist -I/usr/include/scotch -I/usr/include/hdf5/openmpi -I/home/rl420/ICFERST2_new_master/multifluids_icferst/include/ -DHAVE_PETSC  -I/usr/include/vtk-7.1 -I/usr/include/freetype2 -I/usr/include -I/usr/include/x86_64-linux-gnu -I/usr/lib/x86_64-linux-gnu/openmpi/include/openmpi -I/usr/lib/x86_64-linux-gnu/openmpi/include -I/usr/include/python3.8 -I/usr/include/hdf5/openmpi -I/usr/include/jsoncpp -I/usr/include/libxml2 -I/usr/include/tcl   -DHAVE_VTK  -I/usr/include/python3.8 -DHAVE_NUMPY -I/home/rl420/ICFERST1/multifluids_icferst/venv/lib/python3.8/site-packages/numpy/core/include  -O3 -I../include -I/usr/local/include
+
+LINKER  = mpicxx  -Xlinker -export-dynamic -Wl,-O1 -Wl,-Bsymbolic-functions  
+
+MAKE    = make
+
+AR  = ar
+ARFLAGS = cr
+
+OBJS_PREDIAGNOSTICS = Extract_From_State.o multi_matrix_operations.o \
+multi_shape_fct_ND.o multi_shape_fct.o multi_sparsity.o multi_eos.o cv-adv-dif.o multi_solvers.o \
+multi_compositional.o multi_dyncore_dg.o \
+multi_interpolation.o \
+
+OBJS_POSTDIAGNOSTICS = Multiphase_TimeLoop.o Multiphase_Prototype_Wrapper.o
+
+.SUFFIXES: .f90 .F90 .c .cpp .o .a
+
+.f90.o:
+	@echo "        FC $<"
+	@ mkdir -p ../include
+	$(FC) $(FCFLAGS) $(GENFLAGS) -c $<
+.F90.o:
+	@echo "        FC $<"
+	@ mkdir -p ../include
+	$(FC) $(FCFLAGS) $(GENFLAGS) -c $<
+.c.o:
+	@echo "        CC $<"
+	$(CC) $(CFLAGS) -c $<
+.cpp.o:
+	@echo "        CXX $<"
+	$(CXX) $(CXXFLAGS) -c $<
+
+default: $(OBJS_PREDIAGNOSTICS)  $(OBJS_POSTDIAGNOSTICS) 
+
+pre_diagnostics: $(OBJS_PREDIAGNOSTICS)
+
+post_diagnostics: $(OBJS_POSTDIAGNOSTICS)
+
+# cancel implicit rule that tries to make .o from .mod with m2c:
+%.o: %.mod
+
+clean:
+	@rm -f *.o *.d *.mod $(PROG_EXEC)
+
+include Makefile.dependencies

ICFERST/src/multi_eos.F90

@@ -436,6 +436,7 @@ subroutine Calculate_Rho_dRhodP( state, packed_state, iphase, icomp, &
         have_temperature_field = ( stat == 0 )
         Concentration => extract_scalar_field( state( iphase ), 'Concentration', stat )
         have_concentration_field = ( stat == 0 )
+
 
         assert( node_count( pressure ) == size( rho ) )
         assert( node_count( pressure ) == size( drhodp ) )
@@ -587,6 +588,8 @@ subroutine Calculate_Rho_dRhodP( state, packed_state, iphase, icomp, &
               end if
               deallocate( eos_coefs )
 
+
+          !revised version, bypass all concentration calculations when concentration field is not specified!
           elseif( trim( eos_option_path ) == trim( option_path_comp ) // '/BW_eos' ) then
             ! Batzle and Wang (1992) EOS
             ! Use for basin scale simulations.
@@ -596,24 +599,23 @@ subroutine Calculate_Rho_dRhodP( state, packed_state, iphase, icomp, &
             !Calculate the freshwater density
             rho = 1e3 * ( 1 + 1e-6 * (-80*(temperature % val - 273.15) - 3.3*((temperature % val - 273.15))**2 + 0.00175*((temperature % val - 273.15))**3 + 489*pressure%val(1,1,:)*1e-6 - 2*(temperature % val - 273.15)*pressure%val(1,1,:)*1e-6 + 0.016*((temperature % val - 273.15))**2*pressure%val(1,1,:)*1e-6 - 1.3e-5*((temperature % val - 273.15))**3*pressure%val(1,1,:)*1e-6 -0.333*(pressure%val(1,1,:)*1e-6)**2 - 0.002*(temperature % val - 273.15)*(pressure%val(1,1,:)*1e-6)**2))
             !Add the brine contribution
-            rho = rho + 1e3*Concentration % val * (0.668 + 0.44*Concentration % val + 1e-6 * (300*pressure%val(1,1,:)*1e-6 - 2400*pressure%val(1,1,:)*1e-6*Concentration % val + (temperature % val - 273.15) * (80 + 3*(temperature % val - 273.15) - 3300*Concentration % val - 13*pressure%val(1,1,:)*1e-6 + 47*pressure%val(1,1,:)*1e-6*Concentration % val)))
-
+            if (have_concentration_field) rho = rho + 1e3*Concentration % val * (0.668 + 0.44*Concentration % val + 1e-6 * (300*pressure%val(1,1,:)*1e-6 - 2400*pressure%val(1,1,:)*1e-6*Concentration % val + (temperature % val - 273.15) * (80 + 3*(temperature % val - 273.15) - 3300*Concentration % val - 13*pressure%val(1,1,:)*1e-6 + 47*pressure%val(1,1,:)*1e-6*Concentration % val)))
             if (has_boussinesq_aprox) then
               dRhodP = 0.0
             else
               perturbation_pressure = 1.e-5
-
               RhoPlus = 1e3 * ( 1 + 1e-6 * (-80*(temperature % val - 273.15) - 3.3*((temperature % val - 273.15))**2 + 0.00175*((temperature % val - 273.15))**3 + 489*(pressure%val(1,1,:) + perturbation_pressure)*1e-6 - 2*(temperature % val - 273.15)*(pressure%val(1,1,:) + perturbation_pressure)*1e-6 + 0.016*((temperature % val - 273.15))**2*(pressure%val(1,1,:) + perturbation_pressure)*1e-6 - 1.3e-5*((temperature % val - 273.15))**3*(pressure%val(1,1,:) + perturbation_pressure)*1e-6 -0.333*((pressure%val(1,1,:) + perturbation_pressure)*1e-6)**2 - 0.002*(temperature % val - 273.15)*((pressure%val(1,1,:) + perturbation_pressure)*1e-6)**2))
-              !Add the brine contribution
-              RhoPlus = RhoPlus + 1e3*Concentration % val * (0.668 + 0.44*Concentration % val + 1e-6 * (300*(pressure%val(1,1,:) + perturbation_pressure)*1e-6 - 2400*(pressure%val(1,1,:) + perturbation_pressure)*1e-6*Concentration % val + (temperature % val - 273.15) * (80 + 3*(temperature % val - 273.15) - 3300*Concentration % val - 13*(pressure%val(1,1,:) + perturbation_pressure)*1e-6 + 47*(pressure%val(1,1,:) + perturbation_pressure)*1e-6*Concentration % val)))
-
               RhoMinus = 1e3 * ( 1 + 1e-6 * (-80*(temperature % val - 273.15) - 3.3*((temperature % val - 273.15))**2 + 0.00175*((temperature % val - 273.15))**3 + 489*(pressure%val(1,1,:) - perturbation_pressure)*1e-6 - 2*(temperature % val - 273.15)*(pressure%val(1,1,:) - perturbation_pressure)*1e-6 + 0.016*((temperature % val - 273.15))**2*(pressure%val(1,1,:) - perturbation_pressure)*1e-6 - 1.3e-5*((temperature % val - 273.15))**3*(pressure%val(1,1,:) - perturbation_pressure)*1e-6 -0.333*((pressure%val(1,1,:) - perturbation_pressure)*1e-6)**2 - 0.002*(temperature % val - 273.15)*((pressure%val(1,1,:) - perturbation_pressure)*1e-6)**2))
-              !Add the brine contribution
-              RhoMinus = RhoMinus + 1e3*Concentration % val * (0.668 + 0.44*Concentration % val + 1e-6 * (300*(pressure%val(1,1,:) - perturbation_pressure)*1e-6 - 2400*(pressure%val(1,1,:) - perturbation_pressure)*1e-6*Concentration % val + (temperature % val - 273.15) * (80 + 3*(temperature % val - 273.15) - 3300*Concentration % val - 13*(pressure%val(1,1,:) - perturbation_pressure)*1e-6 + 47*(pressure%val(1,1,:) - perturbation_pressure)*1e-6*Concentration % val)))
-
+              if (have_concentration_field) then 
+                  !Add the brine contribution
+                  RhoPlus = RhoPlus + 1e3*Concentration % val * (0.668 + 0.44*Concentration % val + 1e-6 * (300*(pressure%val(1,1,:) + perturbation_pressure)*1e-6 - 2400*(pressure%val(1,1,:) + perturbation_pressure)*1e-6*Concentration % val + (temperature % val - 273.15) * (80 + 3*(temperature % val - 273.15) - 3300*Concentration % val - 13*(pressure%val(1,1,:) + perturbation_pressure)*1e-6 + 47*(pressure%val(1,1,:) + perturbation_pressure)*1e-6*Concentration % val)))
+                  !Add the brine contribution
+                  RhoMinus = RhoMinus + 1e3*Concentration % val * (0.668 + 0.44*Concentration % val + 1e-6 * (300*(pressure%val(1,1,:) - perturbation_pressure)*1e-6 - 2400*(pressure%val(1,1,:) - perturbation_pressure)*1e-6*Concentration % val + (temperature % val - 273.15) * (80 + 3*(temperature % val - 273.15) - 3300*Concentration % val - 13*(pressure%val(1,1,:) - perturbation_pressure)*1e-6 + 47*(pressure%val(1,1,:) - perturbation_pressure)*1e-6*Concentration % val)))
+              end if
               dRhodP = 0.5 * ( RhoPlus - RhoMinus ) / perturbation_pressure
             endif
 
+
         elseif( trim( eos_option_path ) == trim( option_path_python ) ) then
 
             density => extract_scalar_field( state( iphase ), 'Density', stat )
@@ -2627,6 +2629,7 @@ real function retrieve_reference_density(state, packed_state, iphase, icomp, nph
       type (scalar_field) :: sfield
       type (scalar_field), pointer :: pnt_sfield
       type (vector_field), pointer :: position
+
       !Provide input to find out EOS used
       if( icomp > 0 ) then
           eos_option_path = &
@@ -2682,17 +2685,22 @@ real function retrieve_reference_density(state, packed_state, iphase, icomp, nph
       elseif( trim( eos_option_path ) == trim( option_path_comp ) // '/Linear_eos' ) then
         !!$ Den = den0 * ( 1 + alpha * solute mass fraction - beta * DeltaT )
         call get_option( trim( eos_option_path ) // '/reference_density', ref_rho )
+
+
+      !revised version, assign default values for T0, P0 and C0 if these variables are not initiated!
       elseif( trim( eos_option_path ) == trim( option_path_comp ) // '/BW_eos' ) then
         !!$ Reference density is calculated using reference C0, T0, P0.
-        call get_option( trim( eos_option_path ) // '/C0', ref_C0 )
-        call get_option( trim( eos_option_path ) // '/T0', ref_T0 )
-        call get_option( trim( eos_option_path ) // '/P0', ref_P0 )
+        call get_option( trim( eos_option_path ) // '/T0', ref_T0, default=298.)
+        call get_option( trim( eos_option_path ) // '/P0', ref_P0, default=1e5)
+        call get_option( trim( eos_option_path ) // '/C0', ref_C0 , default=0.)
         !Calculate the reference freshwater density
         ref_rho = 1e3 * ( 1 + 1e-6 * (-80*(ref_T0 - 273.15) - 3.3*((ref_T0 - 273.15))**2 + 0.00175*((ref_T0 - 273.15))**3 + 489*ref_P0*1e-6 - 2*(ref_T0 - 273.15)*ref_P0*1e-6 + 0.016*((ref_T0 - 273.15))**2*ref_P0*1e-6 - 1.3e-5*((ref_T0 - 273.15))**3*ref_P0*1e-6 -0.333*(ref_P0*1e-6)**2 - 0.002*(ref_T0 - 273.15)*(ref_P0*1e-6)**2))
-        !Add the reference brine contribution
         ref_rho = ref_rho + 1e3*ref_C0 * (0.668 + 0.44*ref_C0 + 1e-6 * (300*ref_P0*1e-6 - 2400*ref_P0*1e-6*ref_C0 + (ref_T0 - 273.15) * (80 + 3*(ref_T0 - 273.15) - 3300*ref_C0 - 13*ref_P0*1e-6 + 47*ref_P0*1e-6*ref_C0)))
+
+
       else if( trim( eos_option_path ) == trim( option_path_comp ) // '/Temperature_Pressure_correlation' ) then
         call get_option( trim( option_path_comp ) // '/Temperature_Pressure_correlation/rho0', ref_rho)
+
       elseif( trim( eos_option_path ) == trim( option_path_python ) ) then
         !Here the user has to specify by hand under the boussinesq option
         if (.not.have_option(trim( option_path_python ) // '/Boussinesq_approximation/reference_density')) then