diff --git a/FortranCode/AIOMFAC_inout.f90 b/FortranCode/AIOMFAC_inout.f90
index 9ce46f5..3f79406 100644
--- a/FortranCode/AIOMFAC_inout.f90
+++ b/FortranCode/AIOMFAC_inout.f90
@@ -1,8 +1,8 @@
!****************************************************************************************
!* :: Purpose :: *
!* Subroutine providing an interface to AIOMFAC for computing specific liquid phase *
-!* activity coefficients and phase compositions with consideration of a single liquid *
-!* phase as solution. *
+!* activity coefficients and phase compositions with consideration of a single liquid *
+!* phase as solution. *
!* Input is a specific composition point and temperature for a previously initialized *
!* system. Output is listed in outputvars. *
!* *
@@ -27,30 +27,32 @@
!* program. If not, see . *
!* *
!****************************************************************************************
-SUBROUTINE AIOMFAC_inout(inputconc, xinputtype, TKelvin, nspecies, outputvars, outnames, errorflag, warningflag)
+SUBROUTINE AIOMFAC_inout(inputconc, xinputtype, TKelvin, nspecies, outputvars, outputviscvars, outnames, errorflag, warningflag)
USE ModSystemProp
USE ModAIOMFACvar
USE ModCompScaleConversion
-USE ModCalcActCoeff, ONLY : AIOMFAC_calc
+USE ModCalcActCoeff, ONLY : AIOMFAC_calc, DeltaActivities
USE ModSubgroupProp, ONLY : SMWA, SMWC
IMPLICIT NONE
!interface variables:
REAL(8),DIMENSION(nindcomp),INTENT(IN) :: inputconc !inputconc = the concentration of a given input point (e.g., at an experimental data point)
REAL(8),DIMENSION(6,NKNpNGS),INTENT(OUT) :: outputvars !2-D output array with computed compositions and activities for each species; structure is: | mass-frac., mole-frac., molality, act.coeff., activity, ion-indicator | species-no |
+REAL(8),DIMENSION(2),INTENT(OUT) :: outputviscvars !output array for viscosity related values: | viscosity | model sensitivity |
REAL(8),INTENT(IN) :: TKelvin !the input temperature [K]
LOGICAL(4),INTENT(IN) :: xinputtype
INTEGER(4),INTENT(OUT) :: nspecies, errorflag, warningflag
-CHARACTER(LEN=60),DIMENSION(NKNpNGS),INTENT(OUT) :: outnames
+CHARACTER(LEN=*),DIMENSION(NKNpNGS),INTENT(OUT) :: outnames
!--
!local variables:
CHARACTER(LEN=2) :: cn !this assumes a maximum two-digit component number in the system (max. 99); to be adjusted otherwise.
CHARACTER(LEN=3) :: cino
INTEGER(4) :: i, ion_no, ion_indic, nc, NKSinput, NKSinputp1
REAL(8),PARAMETER :: DEPS = 1.1D1*(EPSILON(DEPS))
-REAL(8) :: wtf_cp, xi_cp, mi_cp, actcoeff_cp, a_cp, sum_ms, sum_mions, sum_miMi
+REAL(8) :: wtf_cp, xi_cp, mi_cp, actcoeff_cp, a_cp, sum_ms, sum_miMi, xtolviscosity !, sum_mions
REAL(8),DIMENSION(nelectrol) :: mixingratio, wtfdry
+REAL(8),DIMENSION(nindcomp) :: xinp, dact, dactcoeff, wfrac
!-------------------------------------------------------------------------------------------
!check for debugging:
@@ -64,6 +66,7 @@ SUBROUTINE AIOMFAC_inout(inputconc, xinputtype, TKelvin, nspecies, outputvars, o
errorflag = 0
warningflag = 0
outputvars = 0.0D0
+outputviscvars = 0.0D0
NKSinput = ninput -nneutral
NKSinputp1 = NKSinput+1
wtfdry(1:NKSinput) = 1.0D0
@@ -71,6 +74,12 @@ SUBROUTINE AIOMFAC_inout(inputconc, xinputtype, TKelvin, nspecies, outputvars, o
mixingratio(1:NKSinput) = 1.0D0
mixingratio(NKSinputp1:) = 0.0D0
nspecies = NKNpNGS
+xtolviscosity = 0.0D0
+IF (nneutral < nspecies) THEN
+ calcviscosity = .false. !currently switched off in presence of electrolytes
+ELSE
+ calcviscosity = .true.
+ENDIF
IF (nneutral < 1) THEN !leave the subroutine and indicate a problem to the calling routine
errorflag = 8 !there must be at least one neutral component in the mixture!
@@ -97,6 +106,16 @@ SUBROUTINE AIOMFAC_inout(inputconc, xinputtype, TKelvin, nspecies, outputvars, o
!.....
CALL MassFrac2MoleFracMolality(wtf, XrespSalt, mrespSalt)
CALL AIOMFAC_calc(wtf, TKelvin) !calculate at given mass fraction and temperature
+
+IF (calcviscosity) THEN
+ xinp(1:nindcomp) = XrespSalt(1:nindcomp)
+ CALL DeltaActivities(xinp, TKelvin, dact, dactcoeff)
+ wfrac = wtf
+ wfrac(1) = wfrac(1) + 0.02D0
+ wfrac = wfrac/(1.0D0 + 0.02D0)
+ CALL MassFrac2MoleFracMolality(wfrac, XrespSalt, mrespSalt)
+ xtolviscosity = XrespSalt(1) - xinp(1)
+ENDIF
!.....
!Output of the AIOMFAC calculated values species-wise to array outputvars (ions separately):
@@ -173,22 +192,36 @@ SUBROUTINE AIOMFAC_inout(inputconc, xinputtype, TKelvin, nspecies, outputvars, o
outputvars(6,nc) = REAL(ion_indic, KIND=8) !the indicator if this species is an inorg. ion or not: 0 = not ion, a number > 200 indicates the ion ID from the list
ENDDO
+! viscosity output
+IF (calcviscosity) THEN
+ outputviscvars(1) = LOG10(etamix)
+ outputviscvars(2) = xtolviscosity*deltaetamix*0.4342944819D0 !the factor 0.4342944.. for conversion from deltaetamix which is an ln() value to log_10().
+ELSE
+ outputviscvars(1) = -999.999D0 !negative/unphysical values to signal "property not calculated"
+ outputviscvars(2) = -999.999D0
+ENDIF
+
!check applicable temperature range and state a warning "errorflag" if violated:
!applicable range for electrolyte-containing mixtures (approx.): 288.0 to 309.0 K (298.15 +- 10 K); (strictly valid range would be 298.15 K only)
!applicable range for electrlyte-free mixtures (approx.): 280.0 to 400.0 K
IF (NGS > 0) THEN !electrolyte-containing
IF (TKelvin > 309.0D0 .OR. TKelvin < 288.0D0) THEN !set warning flag
- IF (warningflag == 0) THEN !don't overwrite an existing warning if non-zero
+ IF (warningflag == 0) THEN !do not overwrite an existing warning when non-zero
warningflag = 10
ENDIF
ENDIF
ELSE
IF (TKelvin > 400.0D0 .OR. TKelvin < 280.0D0) THEN !set warning flag
- IF (warningflag == 0) THEN !don't overwrite an existing warning if non-zero
+ IF (warningflag == 0) THEN
warningflag = 11
ENDIF
ENDIF
ENDIF
+IF (.NOT. calcviscosity) THEN
+ IF (warningflag == 0) THEN
+ warningflag = 16
+ ENDIF
+ENDIF
END SUBROUTINE AIOMFAC_inout
! ======================= END =======================================================
\ No newline at end of file
diff --git a/FortranCode/Main_IO_driver.f90 b/FortranCode/Main_IO_driver.f90
index f9b1ea8..0f0833d 100644
--- a/FortranCode/Main_IO_driver.f90
+++ b/FortranCode/Main_IO_driver.f90
@@ -23,7 +23,7 @@
!* Dept. Atmospheric and Oceanic Sciences, McGill University (2013 - present) *
!* *
!* -> created: 2011 (this file) *
-!* -> latest changes: 2019/01/23 *
+!* -> latest changes: 2019/08/21 *
!* *
!* :: License :: *
!* This program is free software: you can redistribute it and/or modify it under the *
@@ -43,293 +43,80 @@ PROGRAM Main_IO_driver
!module variables:
!USE IFPORT !for directory creation
-USE ModSystemProp, ONLY : compname, compsubgroups, compsubgroupsHTML, compsubgroupsTeX, cpname, &
- & errorflagmix, NGS, nindcomp, ninput, NKNpNGS, nneutral, SetSystem, topsubno, waterpresent
-USE ModSubgroupProp, ONLY : SubgroupNames, SubgroupAtoms, subgrname, subgrnameTeX, subgrnameHTML
+USE ModSystemProp, ONLY : errorflagmix, nindcomp, NKNpNGS, SetSystem, topsubno, waterpresent
+USE ModSubgroupProp, ONLY : SubgroupAtoms, SubgroupNames
USE ModMRpart, ONLY : MRdata
USE ModSRunifac, ONLY : SRdata
IMPLICIT NONE
!set preliminary input-related parameters:
INTEGER(4),PARAMETER :: maxpoints = 1000 !limit maximum number of composition points for web-version
-INTEGER(4),PARAMETER :: ninpmax = 50 !set the maximum number of mixture components allowed (preliminary parameter)
+INTEGER(4),PARAMETER :: ninpmax = 51 !set the maximum number of mixture components allowed (preliminary parameter)
!local variables:
-CHARACTER(LEN=2) :: cn !this assumes a maximum two-digit component number in the system (max. 99); to be adjusted otherwise.
-CHARACTER(LEN=4) :: dashes, equalsigns, pluses, VersionNo
+CHARACTER(LEN=4) :: VersionNo
CHARACTER(LEN=20) :: dummy
-CHARACTER(LEN=50) :: subntxt, txtcheck
-CHARACTER(LEN=150) :: horizline, tablehead, tformat, txtn
-CHARACTER(LEN=3000) :: errlogfile, filename, filepath, filepathout, fname, &
- & mixturestring, txtfilein, txtsubs
-CHARACTER(LEN=20),DIMENSION(ninpmax) :: txtarray
-CHARACTER(LEN=60),DIMENSION(ninpmax) :: cpnameinp !list of assigned component names (from input file)
-CHARACTER(LEN=60),DIMENSION(:),ALLOCATABLE :: outnames
-CHARACTER(LEN=400) :: outtxtleft
-INTEGER(4) :: allocstat, cpno, errlogout, errorflagcalc, errorind, i, inpfilesize, &
- & istat, k, nc, ncp, npoints, nspecies, nspecmax, pointi, qty, standardout, &
- & subg, unito, unitx, warningflag, warningind, watercompno
+CHARACTER(LEN=150) :: tformat
+CHARACTER(LEN=3000) :: filename, filepath, filepathout, fname, txtfilein
+CHARACTER(LEN=200),DIMENSION(:),ALLOCATABLE :: cpnameinp !list of assigned component names (from input file)
+CHARACTER(LEN=200),DIMENSION(:),ALLOCATABLE :: outnames
+INTEGER(4) :: allocstat, errorflagcalc, errorind, i, nc, ncp, npoints, nspecies, nspecmax, pointi, &
+ & unito, warningflag, warningind, watercompno
INTEGER(4),DIMENSION(ninpmax) :: px
-INTEGER(4),DIMENSION(ninpmax,topsubno) :: cpsubg !list of input component subgroups and corresponding subgroup quantities
-REAL(8) :: TKelvin, RH
-REAL(8),DIMENSION(maxpoints) :: T_K
-REAL(8),DIMENSION(:),ALLOCATABLE :: inputconc
-REAL(8),DIMENSION(:,:),ALLOCATABLE :: composition, outputvars
+INTEGER(4),DIMENSION(:,:),ALLOCATABLE :: cpsubg !list of input component subgroups and corresponding subgroup quantities
+REAL(8) :: TKelvin
+REAL(8),DIMENSION(:),ALLOCATABLE :: T_K
+REAL(8),DIMENSION(:),ALLOCATABLE :: inputconc, outputviscvars
+REAL(8),DIMENSION(:,:),ALLOCATABLE :: composition, compos2, outputvars, out_viscdata
REAL(8),DIMENSION(:,:,:),ALLOCATABLE :: out_data
-LOGICAL(4) :: ext, fileexists, fileopened, filevalid, verbose, xinputtype
+LOGICAL(4) :: filevalid, verbose, xinputtype
+!--
+!explicit interfaces:
+INTERFACE
+ SUBROUTINE ReadInputFile(filepath, filename, filepathout, ninpmax, maxpoints, unito, verbose, ncp, npoints, &
+ & warningind, errorind, filevalid, cpnameinp, cpsubg, T_K, composition, xinputtype)
+ USE ModSystemProp, ONLY : topsubno
+ CHARACTER(LEN=3000),INTENT(INOUT) :: filepath, filename, filepathout
+ INTEGER(4),INTENT(IN) :: ninpmax, maxpoints
+ INTEGER(4),INTENT(INOUT) :: unito
+ LOGICAL(4),INTENT(IN) :: verbose
+ INTEGER(4),INTENT(OUT) :: ncp, npoints
+ INTEGER(4),INTENT(INOUT) :: warningind, errorind
+ LOGICAL(4),INTENT(OUT) :: filevalid
+ CHARACTER(LEN=200),DIMENSION(:),INTENT(OUT) :: cpnameinp !list of assigned component names (from input file)
+ INTEGER(4),DIMENSION(:,:),INTENT(OUT) :: cpsubg !list of input component subgroups and corresponding subgroup quantities
+ REAL(8),DIMENSION(:),INTENT(OUT) :: T_K !temperature of data points in Kelvin
+ REAL(8),DIMENSION(:,:),INTENT(OUT) :: composition !array of mixture composition points for which calculations should be run
+ LOGICAL(4),INTENT(OUT) :: xinputtype
+ END SUBROUTINE ReadInputFile
+END INTERFACE
!...................................................................................
+
!
!==== INITIALIZATION section =======================================================
!
-VersionNo = "2.21" !AIOMFAC-web version number (change here if minor or major changes require a version number change)
-ncp = 0
+VersionNo = "2.30" !AIOMFAC-web version number (change here if minor or major changes require a version number change)
+verbose = .true. !if true, some debugging information will be printed to the unit "unito" (errorlog file)
nspecmax = 0
-cpsubg = 0
-standardout = 6 !standard output to terminal console window (Fortran standard unit *, 6, or 101)
-errlogout = 71 !output unit for error log file
-cpnameinp = "none" !default component's name
-T_K = 298.15D0 !default temperature in [K]
-dashes = "----"
-equalsigns = "===="
-pluses = "++++"
-unito = errlogout
errorind = 0 !0 means no error found
warningind = 0 !0 means no warnings found
-filevalid = .false.
-fileexists = .false.
-verbose = .true. !if true, some debugging information will be printed to the unit "unito" (errorlog file)
!
!==== INPUT data section ===========================================================
!
!read command line for text-file name (which contains the input parameters to run the AIOMFAC progam):
CALL GET_COMMAND_ARGUMENT(1, txtfilein)
!---
-!txtfilein = './Inputfiles/input_0008.txt' !just use this for debugging with a specific input file, otherwise comment out
+!txtfilein = './Inputfiles/input_1354.txt' !just use this for debugging with a specific input file, otherwise comment out
!---
filepath = ADJUSTL(TRIM(txtfilein))
WRITE(*,*) ""
WRITE(*,*) "MESSAGE from AIOMFAC-web: program started, command line argument 1 = ", filepath
WRITE(*,*) ""
-!extract filepath from the input filepath (which could include a path to the directory):
-k = LEN_TRIM(filepath)
-IF (k < 1) THEN !no file was provided
- WRITE(*,*) ""
- WRITE(*,*) "ERROR from AIOMFAC-web: no input file was provided!"
- WRITE(*,*) "An input file path needs to be provided via command line argument 2"
- WRITE(*,*) ""
- READ(*,*)
- STOP
-ENDIF
-i = INDEX(filepath, "Inputfiles/input")
-IF (i > 0 .AND. i < k) THEN !found a direcory path in UNIX file path style
- filename = filepath(i+11:k)
- filepath = filepath(1:i+10)
-ELSE
- i = INDEX(filepath, "Inputfiles\input")
- IF (i > 0 .AND. i < k) THEN !found a direcory path in WINDOWS file path style
- filename = filepath(i+11:k)
- filepath = filepath(1:i+10)
- ELSE
- filename = filepath(i+11:k)
- filepath = ""
- ENDIF
-ENDIF
-!-----
-!check / create for associated output directory "Outputfiles":
-i = LEN_TRIM(filepath)
-filepathout = TRIM(filepath(1:i-11))//"Outputfiles/"
-!use the name of the input file to create a corresponding output file name:
-i = INDEX(filename, ".txt") !returns starting position of string input within string filename
-!create an error-logfile associated with the input file name:
-errlogfile = "Errorlog_"//filename(i-4:)
-fname = TRIM(filepathout)//TRIM(errlogfile)
-OPEN (NEWUNIT = unito, FILE = fname, STATUS ='UNKNOWN') !unito is the error / logfile unit
-!-----
-!check if file exists and read its content if true:
-fname = TRIM(filepath)//TRIM(filename)
-INQUIRE(FILE = fname, EXIST = fileexists, SIZE = inpfilesize) !inpfilesize is the file size in [bytes]
-!Delete very large files that can only mean uploaded spam content and not actual input:
-IF (fileexists) THEN
- IF (FLOAT(inpfilesize) > 2E6) THEN !cannot be a valid input file
- fname = TRIM(filepath)//TRIM(filename)
- OPEN (NEWUNIT = unitx, FILE = fname, STATUS='OLD')
- CLOSE(unitx, STATUS='DELETE') !close and delete the file
- fileexists = .false.
- ENDIF
-ENDIF
-
-!attempt to read a supposedly valid input file:
-IF (fileexists) THEN
- IF (verbose) THEN
- WRITE(unito,*) ""
- WRITE(unito,*) "MESSAGE from AIOMFAC: input file exists."
- ENDIF
- fname = TRIM(filepath)//TRIM(filename)
- OPEN (NEWUNIT = unitx, FILE = fname, IOSTAT=istat, ACTION='READ', STATUS='OLD')
- IF (istat /= 0) THEN ! an error occurred
- WRITE(unito,*) ""
- WRITE(unito,*) "MESSAGE from AIOMFAC: an error occurred while trying to open the file! IOSTAT: ", istat
- ENDIF
- !validate file as a correct input text-file (no spam):
- READ(unitx,*) dummy, dummy, dummy, txtcheck
- IF (txtcheck(1:11) == "AIOMFAC-web") THEN
- filevalid = .true.
- IF (verbose) THEN
- WRITE(unito,*) ""
- WRITE(unito,*) "MESSAGE from AIOMFAC: input file has passed the first line text validation and will be read."
- ENDIF
- !read input data from file:
- BACKSPACE unitx !jump back to beginning of record (to the beginning of the line)
- READ(unitx,*) dummy, dummy, dummy, dummy, dummy !read first line
- READ(unitx,*) !read empty line 2
- READ(unitx,*) dummy, dummy !read line 3
- READ(unitx,*) dummy !read line 4
- ELSE !file not valid. It will be closed and deleted below
- filevalid = .false.
- IF (verbose) THEN
- WRITE(unito,*) ""
- WRITE(unito,*) "MESSAGE from AIOMFAC: input file does not pass first line text validation and will be deleted."
- ENDIF
- ENDIF
- !loop over mixture components with variable numbers of subgroups to read (using inner loop):
- IF (filevalid) THEN
- IF (verbose) THEN
- WRITE(unito,*) ""
- WRITE(unito,*) "MESSAGE from AIOMFAC: reading component data from input file."
- ENDIF
- DO ncp = 1,ninpmax
- IF (verbose .AND. istat /= 0) THEN
- WRITE(unito,*) ""
- WRITE(unito,*) "MESSAGE from AIOMFAC: file end found in input file at ncp = ", ncp
- ENDIF
- READ(unitx,*,IOSTAT=istat) txtcheck !read only first argument on this line for subsequent check
- IF (txtcheck(1:4) == pluses .OR. istat /= 0) THEN !"++++" indicates end of this components definition part
- EXIT !exit ncp DO-loop
- ELSE !in this case, argument 3 of txtcheck is the component no.:
- BACKSPACE unitx !jump back to beginning of record (to the beginning of the line)
- READ(unitx,*) dummy, dummy, cpno !read line with component's number
- ENDIF
- READ(unitx,*) dummy, dummy, cpnameinp(cpno) !read line with component's name
- IF (LEN_TRIM(cpnameinp(cpno)) < 1) THEN !no component name was assigned, generate a default name
- WRITE(cn,'(I2.2)') cpno
- cpnameinp(cpno) = "cp_"//cn
- ENDIF
- IF (verbose) THEN
- WRITE(unito,*) ""
- WRITE(unito,*) "MESSAGE from AIOMFAC: found a component cpno =", cpno
- ENDIF
- DO !until exit
- READ(unitx,*,IOSTAT=istat) txtcheck !read only first argument on this line
- !check whether another subgroup is present or not
- IF (txtcheck(1:4) == dashes .OR. istat /= 0) THEN !"----" indicates no more subgroups of this component
- EXIT !leave the inner DO-loop
- ELSE !else continue reading the next subgroup
- BACKSPACE unitx
- READ(unitx,*) dummy, dummy, dummy, subg, qty !continue reading line with subgroup no. and corresp. quantity
- cpsubg(cpno,subg) = cpsubg(cpno,subg)+qty
- ENDIF
- ENDDO
- ENDDO
- ncp = ncp-1 !ncp-1 is the number of different components in this mixture
- IF (verbose) THEN
- WRITE(unito,*) ""
- WRITE(unito,*) "MESSAGE from AIOMFAC: component data read."
- WRITE(unito,*) "MESSAGE from AIOMFAC: number of components: ", ncp
- ENDIF
- READ(unitx,*) dummy, dummy !read mixture composition line
- READ(unitx,*) dummy, dummy, i !read mass fraction? line
- IF (i == 1) THEN !composition in mass fractions
- xinputtype = .false.
- ELSE !composition in mole fractions
- xinputtype = .true.
- ENDIF
- READ(unitx,*) dummy, dummy, i !read mole fraction? line
- READ(unitx,*) dummy !read dashes line
- !now read the lines with the mixture composition points:
- !first allocate the composition array with a max of maxpoints points and the number of components present:
- ALLOCATE(composition(maxpoints,ncp), STAT=allocstat)
- composition = 0.0D0 !initialize array
- READ(unitx,*) txtarray(1:ncp) !read header line of composition table
- DO npoints = 1,maxpoints !or until exit
- READ(unitx,*,IOSTAT=istat) txtcheck !read only the first argument on this line
- IF (txtcheck(1:4) == equalsigns .OR. istat /= 0) THEN !"====" indicates no more composition points (and last line of input file)
- EXIT !leave the DO-loop (normal exit point)
- ELSE IF (IACHAR(txtcheck(1:1)) > 47 .AND. IACHAR(txtcheck(1:1)) < 58) THEN !validate that the data is actual intended input and not some sort of text field spam).
- BACKSPACE unitx
- READ(unitx,*) txtcheck, dummy !read only the first two arguments on this line
- IF (IACHAR(dummy(1:1)) > 47 .AND. IACHAR(dummy(1:1)) < 58) THEN !it is a number
- BACKSPACE unitx
- READ(unitx,*) i, T_K(npoints), composition(npoints,2:ncp) !read the temperature in [K] and composition values of the components [2:ncp] into the array
- ELSE
- IF (npoints == 1) THEN
- filevalid = .false. !file is not valid due to incorrect input in text field
- EXIT
- ELSE
- warningind = 31
- EXIT !abnormal exit point of loop due to non-number entries at a certain line in the text field
- ENDIF
- ENDIF
- ELSE
- IF (npoints == 1) THEN
- filevalid = .false. !file is not valid due to incorrect input in text field
- EXIT
- ELSE
- warningind = 31
- EXIT !abnormal exit point of loop due to non-number entries at a certain line in the text field
- ENDIF
- ENDIF
- ENDDO
- npoints = npoints-1 !the number of composition points
- IF (verbose) THEN
- WRITE(unito,*) ""
- WRITE(unito,*) "MESSAGE from AIOMFAC: composition points read."
- WRITE(unito,*) "MESSAGE from AIOMFAC: number of points: ",npoints
- ENDIF
- IF (.NOT. filevalid) THEN
- !close and delete file from server:
- CLOSE(unitx, STATUS = 'DELETE')
- ELSE
- CLOSE(unitx)
- ENDIF
- !assign component 01 its composition values based on the rest of the component's data:
- DO i = 1,npoints
- composition(i,1) = 1.0D0-SUM(composition(i,2:ncp))
- ENDDO
- ENDIF !filevalid
-ELSE
- WRITE(unito,*) ""
- WRITE(unito,*) "ERROR in AIOMFAC: Input file does not exist at expected location: ", TRIM(filename)
- WRITE(unito,*) ""
-ENDIF !fileexists
-
-IF (fileexists) THEN
- IF (.NOT. filevalid) THEN !input file contains errors or is completely invalid (submitted spam etc.)
- IF (verbose) THEN
- WRITE(unito,*) ""
- WRITE(unito,*) "MESSAGE from AIOMFAC: input file did not pass full validation and may be a spam file."
- WRITE(unito,*) "MESSAGE from AIOMFAC: the input file will be deleted to prevent spam files and malicious code from occupying the server."
- WRITE(unito,*) ""
- ENDIF
- errorind = 32
- !close and delete file from server:
- INQUIRE(FILE = fname, OPENED = fileopened)
- IF (fileopened) THEN
- CLOSE(unitx, STATUS = 'DELETE')
- ENDIF
- ELSE
- !check for any valid points for calculations before initializing AIOMFAC:
- IF (npoints < 1) THEN !no points input for calculations
- errorind = 33 !no input in text field..
- IF (verbose) THEN
- WRITE(unito,*) ""
- WRITE(unito,*) "MESSAGE from AIOMFAC: no composition points have been entered in the text field. There is nothing to calculate."
- WRITE(unito,*) ""
- ENDIF
- filevalid = .false.
- ENDIF
- ENDIF
-ENDIF
-
-IF (fileexists .AND. filevalid) THEN
+ALLOCATE(cpsubg(ninpmax,topsubno), cpnameinp(ninpmax), composition(maxpoints,ninpmax), T_K(maxpoints), STAT=allocstat)
+!--
+CALL ReadInputFile(filepath, filename, filepathout, ninpmax, maxpoints, unito, verbose, ncp, npoints, &
+ & warningind, errorind, filevalid, cpnameinp, cpsubg, T_K, composition, xinputtype)
+!--
+IF (filevalid) THEN
!
!==== AIOMFAC initialization and calculation section ===============================
!
@@ -352,6 +139,12 @@ PROGRAM Main_IO_driver
IF (waterpresent) THEN
watercompno = MAXLOC(cpsubg(1:ncp,16), DIM=1) !usually = 1
ENDIF
+ !transfer composition data to adequate array size:
+ ALLOCATE(compos2(npoints,ncp), STAT=allocstat)
+ DO nc = 1,ncp
+ compos2(1:npoints,nc) = composition(1:npoints,nc)
+ ENDDO
+ DEALLOCATE(cpsubg, composition, STAT=allocstat)
IF (errorflagmix /= 0) THEN !a mixture-related error occured:
CALL RepErrorWarning(unito, errorflagmix, warningflag, errorflagcalc, i, errorind, warningind)
@@ -359,9 +152,10 @@ PROGRAM Main_IO_driver
IF (errorind == 0) THEN !perform AIOMFAC calculations; else jump to termination section
!--
- ALLOCATE(inputconc(nindcomp), outputvars(6,NKNpNGS), outnames(NKNpNGS), out_data(7,npoints,NKNpNGS), STAT=allocstat)
+ ALLOCATE(inputconc(nindcomp), outputvars(6,NKNpNGS), outputviscvars(nindcomp), outnames(NKNpNGS), out_data(7,npoints,NKNpNGS), out_viscdata(3,npoints), STAT=allocstat)
inputconc = 0.0D0
out_data = 0.0D0
+ out_viscdata = 0.0D0
!--
IF (verbose) THEN
WRITE(unito,*) ""
@@ -371,10 +165,10 @@ PROGRAM Main_IO_driver
px = 0
!set AIOMFAC input and call the main AIOMFAC subroutines for all composition points:
DO pointi = 1,npoints !loop over points, changing composition / temperature
- inputconc(1:ncp) = composition(pointi,1:ncp)
+ inputconc(1:ncp) = compos2(pointi,1:ncp)
TKelvin = T_K(pointi)
!--
- CALL AIOMFAC_inout(inputconc, xinputtype, TKelvin, nspecies, outputvars, outnames, errorflagcalc, warningflag)
+ CALL AIOMFAC_inout(inputconc, xinputtype, TKelvin, nspecies, outputvars, outputviscvars, outnames, errorflagcalc, warningflag)
!--
IF (warningflag > 0 .OR. errorflagcalc > 0) THEN
!$OMP CRITICAL errwriting
@@ -386,13 +180,20 @@ PROGRAM Main_IO_driver
DO nc = 1,nspecmax !loop over species (ions dissociated and treated as individual species):
out_data(1:6,pointi,nc) = outputvars(1:6,nc) !out_data general structure: | data columns 1:7 | data point | component no.|
out_data(7,pointi,nc) = REAL(errorflagcalc, KIND=8)
+ out_viscdata(3,pointi) = REAL(errorflagcalc, KIND=8)
IF (errorflagcalc == 0 .AND. warningflag > 0) THEN !do not overwrite an errorflag if present!
- out_data(7,pointi,nc) = REAL(warningflag, KIND=8)
+ IF (warningflag == 16) THEN !a warning that only affects viscosity calc.
+ out_viscdata(3,pointi) = REAL(warningflag, KIND=8)
+ ELSE
+ out_data(7,pointi,nc) = REAL(warningflag, KIND=8)
+ out_viscdata(3,pointi) = REAL(warningflag, KIND=8)
+ ENDIF
ENDIF
IF (px(nc) == 0 .AND. out_data(6,pointi,nc) >= 0.0D0) THEN
px(nc) = pointi !use a point for which this component's abundance is given, i.e. mole fraction(nc) > 0.0!
ENDIF
ENDDO !nc
+ out_viscdata(1:2,pointi) = outputviscvars(1:2) !out_viscdata general structure: | data columns 1:2 | data point |
ENDDO !pointi
!
!==== OUTPUT data-to-file section ==================================================
@@ -410,17 +211,22 @@ PROGRAM Main_IO_driver
WRITE(unito,'(A80)') "................................................................................"
!create an output ASCII text file with an overall mixture header and individual tables for all components / species (in case of ions)
fname = TRIM(filepathout)//TRIM(filename)
- CALL OutputTXT(fname, VersionNo, cpnameinp(1:nspecmax), nspecmax, npoints, watercompno, T_K(1:npoints), px(1:nspecmax), out_data)
+ CALL OutputTXT(fname, VersionNo, cpnameinp(1:nspecmax), nspecmax, npoints, watercompno, T_K(1:npoints), px(1:nspecmax), out_data, out_viscdata)
+ !--
+ !>> write output HTML-file
+ i = LEN_TRIM(filename)
+ filename = filename(1:i-3)//"html"
+ fname = TRIM(filepathout)//TRIM(filename)
+ CALL OutputHTML(fname, VersionNo, cpnameinp(1:nspecmax), nspecmax, npoints, watercompno, T_K(1:npoints), px(1:nspecmax), out_data, out_viscdata)
!
!==== TERMINATION section ==========================================================
!
- DEALLOCATE(inputconc, outputvars, outnames, out_data, STAT=allocstat)
- ext = ALLOCATED(composition)
- IF (ext) THEN
- DEALLOCATE(composition, STAT=allocstat)
+ DEALLOCATE(inputconc, outputvars, outputviscvars, outnames, out_data, out_viscdata, T_K, cpnameinp, STAT=allocstat)
+ IF (ALLOCATED(compos2)) THEN
+ DEALLOCATE(compos2, STAT=allocstat)
ENDIF
ENDIF !errorind
-ENDIF !fileexists and valid
+ENDIF !file valid
WRITE(unito,*) "+-+-+-+-+"
WRITE(unito,*) "Final warning indicator (an entry '00' means no warnings found):"
@@ -440,4 +246,4 @@ PROGRAM Main_IO_driver
!
!==== THE END ======================================================================
!
-END PROGRAM Main_IO_driver
+END PROGRAM Main_IO_driver
\ No newline at end of file
diff --git a/FortranCode/ModAIOMFACvar.f90 b/FortranCode/ModAIOMFACvar.f90
index a1f68e6..6fbdc15 100644
--- a/FortranCode/ModAIOMFACvar.f90
+++ b/FortranCode/ModAIOMFACvar.f90
@@ -11,8 +11,8 @@
!* *
!* -> created: 2009 (based on non-module version from 2005) *
!* -> latest changes: 2018/05/29 *
-!* *
-!* :: License :: *
+!* *
+!* :: License :: *
!* This program is free software: you can redistribute it and/or modify it under the *
!* terms of the GNU General Public License as published by the Free Software *
!* Foundation, either version 3 of the License, or (at your option) any later *
@@ -33,22 +33,22 @@ MODULE ModAIOMFACvar
IMPLICIT NONE
!Public module variables:
-REAL(8),PUBLIC :: alphaHSO4, diffKHSO4, ionicstrength, lastTK, meanSolventMW, &
+REAL(8),PUBLIC :: alphaHSO4, deltaetamix, diffKHSO4, etamix, ionicstrength, lastTK, meanSolventMW, &
& SumIonMolalities, T_K, Tmolal, TmolalSolvMix, Xwdissoc
-REAL(8),DIMENSION(:),ALLOCATABLE,PUBLIC :: actcoeff_a, actcoeff_c, actcoeff_n, activity, galrln, gamrln, &
+REAL(8),DIMENSION(:),ALLOCATABLE,PUBLIC :: actcoeff_a, actcoeff_c, actcoeff_n, activity, eta0, eta_cpn, fragil, galrln, gamrln, &
& gasrln, gclrln, gcmrln, gcsrln, gnlrln, gnmrln, gnsrln, ionactivityprod, lnactcoeff_a, &
& lnactcoeff_c, lnactcoeff_n, lnmeanmactcoeff, meanmolalactcoeff, mrespSalt, SMA, SMC, solvmixcorrMRa, &
- & solvmixcorrMRc, wtf, X, XN, XrespSalt
+ & solvmixcorrMRc, Tglass0, wtf, X, XN, XrespSalt
REAL(8),DIMENSION(201:261),PUBLIC :: actcoeff_ion, molality_ion
LOGICAL(4),PUBLIC :: DebyeHrefresh
!..................................................
!make all variables of this module threadprivate for use in parallel execution with openMP:
-!$OMP THREADPRIVATE( alphaHSO4, diffKHSO4, lastTK, meanSolventMW, SumIonMolalities, &
- !$OMP & T_K, Xwdissoc, actcoeff_a, actcoeff_c, actcoeff_n, activity, galrln, gamrln, gasrln, &
+!$OMP THREADPRIVATE( alphaHSO4, deltaetamix, diffKHSO4, etamix, lastTK, meanSolventMW, SumIonMolalities, &
+ !$OMP & T_K, Xwdissoc, actcoeff_a, actcoeff_c, actcoeff_n, activity, eta0, eta_cpn, fragil, galrln, gamrln, gasrln, &
!$OMP & gclrln, gcmrln, gcsrln, gnlrln, gnmrln, gnsrln, ionactivityprod, ionicstrength, lnactcoeff_a, &
!$OMP & lnactcoeff_c, lnactcoeff_n, lnmeanmactcoeff, meanmolalactcoeff, mrespSalt, SMA, SMC, &
- !$OMP & solvmixcorrMRa, solvmixcorrMRc, Tmolal, TmolalSolvMix, wtf, X, XN, XrespSalt, DebyeHrefresh )
+ !$OMP & solvmixcorrMRa, solvmixcorrMRc, Tglass0, Tmolal, TmolalSolvMix, wtf, X, XN, XrespSalt, DebyeHrefresh )
!==========================================================================================================================
CONTAINS
@@ -64,13 +64,14 @@ SUBROUTINE AllocModAIOMFACvar()
!-- allocate several composition-dependent variables:
IF (ALLOCATED(wtf)) THEN
DEALLOCATE ( sma, smc, wtf, X, XN, solvmixcorrMRc, solvmixcorrMRa, XrespSalt, mrespSalt, &
- & activity, meanmolalactcoeff, actcoeff_n, actcoeff_c, actcoeff_a, ionactivityprod, &
+ & activity, meanmolalactcoeff, actcoeff_n, actcoeff_c, actcoeff_a, ionactivityprod, eta0, eta_cpn, fragil, &
& gnlrln, gclrln, galrln, gnmrln, gcmrln, gamrln, gnsrln, gcsrln, gasrln, lnactcoeff_n, lnactcoeff_c, &
- & lnactcoeff_a, lnmeanmactcoeff )
+ & lnactcoeff_a, lnmeanmactcoeff, Tglass0 )
ENDIF
ALLOCATE( sma(NGI), smc(NGI), wtf(nindcomp), X(NKNpNGS), XN(NKNpNGS), XrespSalt(nindcomp), &
& mrespSalt(nindcomp), activity(nindcomp), meanmolalactcoeff(nelectrol), actcoeff_n(nneutral), actcoeff_c(NGI), &
- & actcoeff_a(NGI), ionactivityprod(nelectrol), solvmixcorrMRc(NGI), solvmixcorrMRa(NGI), &
+ & actcoeff_a(NGI), ionactivityprod(nelectrol), solvmixcorrMRc(NGI), solvmixcorrMRa(NGI), eta0(NKNpNGS), fragil(nindcomp), &
+ & eta_cpn(nindcomp), Tglass0(nindcomp), &
& lnactcoeff_n(nneutral), gnlrln(nneutral), gnmrln(nneutral), gnsrln(nneutral), &
& lnactcoeff_c(NGI), gclrln(NGI), gcmrln(NGI), gcsrln(NGI), lnmeanmactcoeff(nelectrol), &
& lnactcoeff_a(NGI), galrln(NGI), gamrln(NGI), gasrln(NGI) )
diff --git a/FortranCode/ModCalcActCoeff.f90 b/FortranCode/ModCalcActCoeff.f90
index 2dd86af..77a2199 100644
--- a/FortranCode/ModCalcActCoeff.f90
+++ b/FortranCode/ModCalcActCoeff.f90
@@ -11,8 +11,8 @@
!* *
!* -> created: 2018 *
!* -> latest changes: 2018/05/28 *
-!* *
-!* :: License :: *
+!* *
+!* :: License :: *
!* This program is free software: you can redistribute it and/or modify it under the *
!* terms of the GNU General Public License as published by the Free Software *
!* Foundation, either version 3 of the License, or (at your option) any later *
@@ -31,6 +31,8 @@
!* - SUBROUTINE HSO4DissociationCheck *
!* CONTAINS: - SUBROUTINE DiffKsulfuricDissoc *
!* - FUNCTION fHSO4dissoc *
+!* - SUBROUTINE DeltaActivities *
+!* - SUBROUTINE partialdactcoeff
!* *
!****************************************************************************************
MODULE ModCalcActCoeff
@@ -115,7 +117,7 @@ SUBROUTINE AIOMFAC_calc(WTFin, TKelvin)
actcoeff_a = 0.0D0
actcoeff_ion = 0.0D0 !to save activity coefficients of ions at their ID number
molality_ion = 0.0D0
-
+
!ln of the activity coefficient for the neutrals:
DO CONCURRENT (I = 1:nneutral) !DO I = 1,nneutral
IF (wtf(I) > dtiny) THEN
@@ -651,5 +653,205 @@ END FUNCTION fHSO4dissoc
END SUBROUTINE HSO4DissociationCheck
!==========================================================================================================================
-
+
+
+ !****************************************************************************************
+ !* :: Purpose :: *
+ !* Subroutine to calculate total derivatives (finite differences) of the AIOMFAC *
+ !* activities and activity coefficients at a given composition (xin) of mixture nd. *
+ !* *
+ !* :: Author & Copyright :: *
+ !* Andi Zuend, *
+ !* Dept. Chem. Engineering, California Institute of Technology (2009 - 2012), *
+ !* Dept. Atmospheric and Oceanic Sciences, McGill University *
+ !* *
+ !* -> created: 2010 *
+ !* -> latest changes: 2018/05/31 *
+ !* *
+ !****************************************************************************************
+ SUBROUTINE DeltaActivities(xin, TKelvin, dact, dactcoeff)
+
+ USE ModSystemProp, ONLY : nindcomp
+ USE ModAIOMFACvar, ONLY : deltaetamix
+
+ IMPLICIT NONE
+
+ !interface variables:
+ REAL(8),DIMENSION(nindcomp),INTENT(IN) :: xin
+ REAL(8),DIMENSION(nindcomp),INTENT(OUT) :: dact, dactcoeff
+ REAL(8),INTENT(IN) :: TKelvin
+ !local variables:
+ INTEGER(4) :: ni, j
+ LOGICAL(4) :: xinputtype
+ REAL(8),DIMENSION(nindcomp) :: xinit
+ REAL(8),DIMENSION(nindcomp,nindcomp) :: partdact, partdactcoeff
+ REAL(8) :: Dtiny, dn, partdact_ji, partdactcoeff_ji
+ PARAMETER (Dtiny = 1.79D1*EPSILON(1.0D0))
+ PARAMETER (dn = 1.0D-9) ![mol] a small, but not too small molar (change) for the numerical differentiation
+ PARAMETER (xinputtype = .true.) !mole fraction input concentration
+ !...........................................................
+ !parameters:
+ xinit = xin !xinit is locally stored and unaffected by changes when reentrant code will be called that might feed back on xin!
+
+ !calculate the finite differences (numerical "partial derivatives") of the component ni's activity at composition xinit while holding the moles of the other componentes fixed.
+ !compositions for the forward / backward differences at the given point xinit (component nnvar is the one to be enhanced/diminished):
+ DO ni = 1,nindcomp !loop
+ !calling the subroutine calculating the partial derivative of the activity and act. coeff. of component j with respect to moles of component ni:
+ DO j = 1,nindcomp
+ CALL partialdactcoeff(xinit, TKelvin, j, ni, partdact_ji, partdactcoeff_ji)
+ partdact(j,ni) = partdact_ji
+ partdactcoeff(j,ni) = partdactcoeff_ji
+ ENDDO
+ ENDDO !ni
+ !---
+ !calculate the total derivative for each component from the absolute values of the partial derivatives:
+ dact(1:nindcomp) = 0.0D0
+ dactcoeff(1:nindcomp) = 0.0D0
+ DO ni = 1,nindcomp
+ dact(ni) = SUM(ABS(partdact(ni,1:nindcomp))) !we skip here the *dn (which would allow calculating the total differential, as it leads to small numbers and sensitivity will be normalized later anyways)
+ dactcoeff(ni) = SUM(ABS(partdactcoeff(ni,1:nindcomp)))
+ ENDDO
+ !for viscosity sensitivity estimate as a function of a change in mole fraction of component 1 (water):
+ deltaetamix = ABS(deltaetamix)
+
+ END SUBROUTINE DeltaActivities
+ !==========================================================================================================================
+
+
+ !****************************************************************************************
+ !* :: Purpose :: *
+ !* Subroutine to calculate partial derivatives (finite differences) of the AIOMFAC *
+ !* activities and activity coefficients of a certain mixture component "j" with *
+ !* respect to a tiny change in moles of component "i" at a composition vector (xin). *
+ !* *
+ !* :: Author & Copyright :: *
+ !* Andi Zuend, *
+ !* Dept. Chem. Engineering, California Institute of Technology (2009 - 2012), *
+ !* Dept. Atmospheric and Oceanic Sciences, McGill University *
+ !* *
+ !* -> created: 2010 *
+ !* -> latest changes: 2018/05/31 *
+ !* *
+ !****************************************************************************************
+ SUBROUTINE partialdactcoeff(xin, TKelvin, j, i, partdact_ji, partdactcoeff_ji)
+
+ USE ModSystemProp, ONLY : nindcomp, nneutral, nd, Mmass, errorflagcalc
+ USE ModCompScaleConversion
+ USE ModAIOMFACvar, ONLY : activity, meanmolalactcoeff, actcoeff_n, deltaetamix, etamix
+
+ IMPLICIT NONE
+
+ !interface variables:
+ INTEGER(4),INTENT(IN) :: j, i
+ REAL(8),DIMENSION(nindcomp),INTENT(IN) :: xin
+ REAL(8),INTENT(IN) :: TKelvin
+ REAL(8),INTENT(OUT) :: partdact_ji, partdactcoeff_ji
+ !local variables:
+ LOGICAL(4) :: xinputtype
+ REAL(8),DIMENSION(nindcomp) :: xinit, xplus, wtf
+ REAL(8) :: ntplus, Dtiny, dn, actplus, actinit, actcoeffplus, actcoeffinit, etamixplus, etamixinit
+ PARAMETER (Dtiny = 1.11D0*EPSILON(1.0D0))
+ PARAMETER (dn = 1.0D-9) ![mol] a small molar (change) for the numerical differentiation
+ PARAMETER (xinputtype = .true.) !mole fraction input concentration
+ !...........................................................
+ !parameters:
+ xinit = xin !xinit is locally stored and unaffected by changes when reentrant code will be called that might feed back on xin!
+
+ !calculate the finite differences (numerical "partial derivatives") of component j's activity and activity coeff. with respect to component i,
+ !at composition xinit, while holding the moles of the other componentes fixed.
+ !compositions for the forward / backward differences at the given point xinit (component i is the one to be enhanced/diminished):
+ IF (j <= nindcomp .AND. i <= nindcomp) THEN !the nindcomp comparisons are necessary for the cases where HSO4- dissociation can change the nindcomp during calc.
+ ntplus = 1.0D0+dn !as the molar sum of all components can be set arbitrarily to 1.0
+ !calculate the mole fractions at the slightly changed compositions:
+ xplus(1:i-1) = xinit(1:i-1)/ntplus
+ xplus(i) = (xinit(i)+dn)/ntplus
+ xplus(i+1:nindcomp) = xinit(i+1:nindcomp)/ntplus
+ !..
+ !Call AIOMFAC_calc to calculate the activity at composition xplus:
+ CALL MoleFrac2MassFrac(xplus, Mmass, wtf)
+ CALL AIOMFAC_calc(wtf, TKelvin)
+ actplus = activity(j)
+ IF (j <= nneutral) THEN !neutral component
+ actcoeffplus = actcoeff_n(j)
+ ELSE !electrolyte component
+ actcoeffplus = meanmolalactcoeff(j-nneutral)
+ ENDIF
+ IF (j == 1 .AND. i == 1) THEN !save the etamix value for the sensitivity of viscosity
+ etamixplus = etamix
+ ENDIF
+ !..
+ !calculate the activity & activity coeff. at the initial composition xinit:
+ CALL MoleFrac2MassFrac(xinit, Mmass, wtf)
+ CALL AIOMFAC_calc(wtf, TKelvin)
+ actinit = activity(j)
+ IF (j <= nneutral) THEN
+ actcoeffinit = actcoeff_n(j)
+ ELSE
+ actcoeffinit = meanmolalactcoeff(j-nneutral)
+ ENDIF
+ IF (j == 1 .AND. i == 1) THEN !save the etamix value for the sensitivity of viscosity
+ etamixinit = etamix
+ ENDIF
+ !..
+ !calculate the numerical forward differences with respect to activity at point xinit
+ !(partial derivative of component j with resp. to a small molar change, dn, in component i)
+ IF (actplus > 0.0D0 .AND. actinit > 0.0D0) THEN
+ IF (ABS(actplus-actinit) < 1.0D-292) THEN !floating underflow risk
+ IF (actplus-actinit < 0.0D0) THEN !"negative zero"
+ partdact_ji = -1.0D-292 !almost zero
+ ELSE !"positive zero"
+ partdact_ji = 1.0D-292 !almost zero
+ ENDIF
+ ELSE !no underflow
+ IF (actplus -actinit > 1.0D292) THEN !floating overflow risk
+ partdact_ji = 1.0D292+(LOG(actplus)-LOG(actinit)) !set to huge positive number to allow following computations but prevent overflow
+ ELSE IF (actplus -actinit < -1.0D292) THEN !floating overflow risk
+ partdact_ji = -1.0D292-(LOG(actplus)-LOG(actinit)) !set to huge negative number to allow following computations but prevent overflow
+ ELSE
+ partdact_ji = (actplus -actinit)/dn !the numerical derivatives with respect to component i (while all other component's moles are kept const.)
+ ENDIF
+ ENDIF
+ IF (ABS(actcoeffplus -actcoeffinit) < 1.0D-292) THEN !floating underflow risk
+ IF (actcoeffplus -actcoeffinit < 0.0D0) THEN !"negative zero"
+ partdactcoeff_ji = -1.0D-292 !almost zero
+ ELSE !"positive zero"
+ partdactcoeff_ji = 1.0D-292 !almost zero
+ ENDIF
+ ELSE !ok, no underflow
+ IF (actcoeffplus -actcoeffinit > 1.0D292) THEN !floating overflow risk
+ partdactcoeff_ji = 1.0D292+(LOG(actcoeffplus)-LOG(actcoeffinit)) !set to huge positive number to allow following computations but prevent overflow
+ ELSE IF (actcoeffplus -actcoeffinit < -1.0D292) THEN !floating overflow risk
+ partdactcoeff_ji = -1.0D292-(LOG(actcoeffplus)-LOG(actcoeffinit)) !set to huge negative number to allow following computations but prevent overflow
+ ELSE
+ partdactcoeff_ji = (actcoeffplus -actcoeffinit)/dn
+ ENDIF
+ ENDIF
+ ELSE IF (xinit(j) > Dtiny .AND. (actplus < Dtiny .OR. actinit < Dtiny)) THEN !there is some amount of a component in the mixture, but the model activity coefficient is too large (thus, very steep derivative)
+ partdact_ji = 1.11111111D5 !a specific value indicating a floating point overflow in AIOMFAC_calc, while not suppressing the feedback of such a value.
+ partdactcoeff_ji = 1.11111111D5
+ errorflagcalc = 2
+ ELSE !exceptions where problem occurred
+ partdact_ji = 0.0D0
+ partdactcoeff_ji = 0.0D0
+ errorflagcalc = 3
+ ENDIF
+ !**--
+ IF (j == 1 .AND. i == 1) THEN
+ !partial forward difference for mixture viscosity:
+ IF (etamixinit > 0.0D0) THEN
+ deltaetamix = (LOG(etamixplus) -LOG(etamixinit))/dn !the numerical derivatives with respect to component i (while all other component's moles are kept const.)
+ ELSE
+ deltaetamix = 0.0D0
+ ENDIF
+ ENDIF
+ !**--
+ ELSE
+ partdact_ji = 0.0D0
+ partdactcoeff_ji = 0.0D0
+ errorflagcalc = 4
+ ENDIF
+
+ END SUBROUTINE partialdactcoeff
+ !==========================================================================================================================
+
END MODULE ModCalcActCoeff
\ No newline at end of file
diff --git a/FortranCode/ModComponentNames.f90 b/FortranCode/ModComponentNames.f90
index b5afd12..4b61d81 100644
--- a/FortranCode/ModComponentNames.f90
+++ b/FortranCode/ModComponentNames.f90
@@ -11,8 +11,8 @@
!* *
!* -> created: 2006 *
!* -> latest changes: 2018/05/22 *
-!* *
-!* :: License :: *
+!* *
+!* :: License :: *
!* This program is free software: you can redistribute it and/or modify it under the *
!* terms of the GNU General Public License as published by the Free Software *
!* Foundation, either version 3 of the License, or (at your option) any later *
@@ -776,8 +776,8 @@ PURE SUBROUTINE names_mix(ncomp, nneutral, CompN, compname, compnameTeX, ionname
INTEGER(4),INTENT(IN) :: ncomp, nneutral !the component number (internal numbering as the compN in definemixtures, LRdata, etc.)
INTEGER(4),DIMENSION(ncomp),INTENT(IN) :: CompN
REAL(8),DIMENSION(ncomp),INTENT(OUT) :: OtoCratio, HtoCratio
- CHARACTER(LEN=60),DIMENSION(ncomp),INTENT(OUT) :: compname, compnameTeX
- CHARACTER(LEN=16),DIMENSION(NGS),INTENT(OUT) :: ionname, ionnameTeX
+ CHARACTER(LEN=*),DIMENSION(ncomp),INTENT(OUT) :: compname, compnameTeX
+ CHARACTER(LEN=*),DIMENSION(NGS),INTENT(OUT) :: ionname, ionnameTeX
!local vars:
CHARACTER(LEN=16) :: txt
INTEGER(4) :: i, k, cnt, cn, an
diff --git a/FortranCode/ModMRpart.f90 b/FortranCode/ModMRpart.f90
index 81cd4b0..c54ebe1 100644
--- a/FortranCode/ModMRpart.f90
+++ b/FortranCode/ModMRpart.f90
@@ -12,9 +12,9 @@
!* Dept. Atmospheric and Oceanic Sciences, McGill University *
!* *
!* -> created: 2004 *
-!* -> latest changes: 2018/05/26 *
-!* *
-!* :: License :: *
+!* -> latest changes: 2018/08/22 *
+!* *
+!* :: License :: *
!* This program is free software: you can redistribute it and/or modify it under the *
!* terms of the GNU General Public License as published by the Free Software *
!* Foundation, either version 3 of the License, or (at your option) any later *
@@ -87,6 +87,8 @@ MODULE ModMRpart
DATA cTABna(1:75,02) / 5.888996D-02, 1.363867D-01, 3.860389D-02, -7.777778D+04, -2.249250D-02, -7.777778D+04, 0.000000D+00, 3.767145D-02, 1.363866D-01, 5.886737D-02, 1.362981D-01, -7.777778D+04, 4.161487D-02, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, 3.009678D-02, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, 3.009678D-02, -2.249250D-02, -7.777778D+04, 3.009678D-02, -2.249250D-02, 3.009678D-02, 5.888996D-02, 5.888996D-02, 3.656014D-02, -2.249250D-02, -7.777778D+04, -7.777778D+04, 1.912237D-02, 3.009678D-02, 8.322975D-02, -7.777778D+04 /
DATA bTABna(1:75,03) / 4.136791D-02, 8.229121D-02, 4.128208D-02, -7.777778D+04, 4.686682D-03, -7.777778D+04, 0.000000D+00, 2.998976D-02, 6.480397D-02, 3.717467D-02, 6.475417D-02, -7.777778D+04, 4.133287D-02, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, 2.998974D-02, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, 2.998974D-02, 4.686682D-03, -7.777778D+04, 2.998974D-02, 4.686682D-03, 2.998974D-02, 4.136791D-02, 4.136791D-02, 4.094891D-02, 4.686682D-03, -7.777778D+04, -7.777778D+04, 4.601955D-02, 2.998974D-02, 8.266574D-02, -7.777778D+04 /
DATA cTABna(1:75,03) / 3.439384D-02, 5.607449D-02, 3.401668D-02, -7.777778D+04, 1.476176D-03, -7.777778D+04, 0.000000D+00, 1.237431D-02, 4.151743D-02, 1.807473D-02, 4.150937D-02, -7.777778D+04, 2.391176D-02, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, 1.237430D-02, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, 1.237430D-02, 1.476176D-03, -7.777778D+04, 1.237430D-02, 1.476176D-03, 1.237430D-02, 3.439384D-02, 3.439384D-02, 1.688331D-02, 1.476176D-03, -7.777778D+04, -7.777778D+04, 2.538794D-02, 1.237430D-02, 4.782352D-02, -7.777778D+04 /
+DATA bTABna(1:75,04) / -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, 0.000000D+00, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04 /
+DATA cTABna(1:75,04) / -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, 0.000000D+00, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04 /
DATA bTABna(1:75,05) / 4.233234D-02, 6.701021D-02, 4.209408D-02, -7.777778D+04, -2.707269D-02, -7.777778D+04, 0.000000D+00, 4.535367D-03, 4.542737D-02, 4.139841D-02, 4.542735D-02, -7.777778D+04, 4.773836D-03, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -2.755684D-03, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -2.755684D-03, -2.707269D-02, -7.777778D+04, -2.755684D-03, -2.707269D-02, -2.755684D-03, 4.233234D-02, 4.233234D-02, 3.000031D-02, -2.707269D-02, -7.777778D+04, -7.777778D+04, -2.229886D-02, -2.755684D-03, 9.547671D-03, -7.777778D+04 /
DATA cTABna(1:75,05) / 4.853253D-02, 1.230290D-01, 4.848432D-02, -7.777778D+04, 8.326951D-03, -7.777778D+04, 0.000000D+00, 2.880167D-02, 6.624150D-02, 4.875210D-02, 6.624148D-02, -7.777778D+04, 8.524851D-02, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, 6.663517D-04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, 6.663517D-04, 8.326951D-03, -7.777778D+04, 6.663517D-04, 8.326951D-03, 6.663517D-04, 4.853253D-02, 4.853253D-02, 3.967023D-02, 8.326951D-03, -7.777778D+04, -7.777778D+04, 9.357546D-02, 6.663517D-04, 1.704970D-01, -7.777778D+04 /
DATA bTABna(1:75,08) / 1.244526D-01, 2.299447D-01, 9.637726D-02, -7.777778D+04, -6.100228D-02, -7.777778D+04, 0.000000D+00, -7.777778D+04, 2.333990D-01, 1.089460D-01, 2.333990D-01, -7.777778D+04, 6.071562D-02, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, 7.999459D-02, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, -7.777778D+04, 7.999459D-02, -6.100228D-02, -7.777778D+04, 7.999459D-02, -6.100228D-02, 7.999459D-02, 1.244526D-01, 1.244526D-01, 9.617478D-02, -6.100228D-02, -7.777778D+04, -7.777778D+04, -2.866620D-04, 7.999459D-02, 1.214312D-01, -7.777778D+04 /
@@ -224,13 +226,13 @@ SUBROUTINE LR_MR_activity()
lrw1 = -0.5D0/SI2
DO CONCURRENT (J = 1:NGI)
oexp1(1:NGI) = omega(1:NGI,J)*SI2
- WHERE (oexp1(1:NGI) > 300.0D0) !prevent floating point underflow as the second term becomes practically zero:
+ WHERE (oexp1(1:NGI) > 300.0D0) !prevent floating point underflow as the second term becomes practically zero:
Bca(1:NGI,J) = bac(1:NGI,J)
ELSEWHERE
Bca(1:NGI,J) = bac(1:NGI,J) +cac(1:NGI,J)*EXP(-oexp1(1:NGI))
ENDWHERE
oexp2(1:NGI) = omega2(1:NGI,J)*SI2
- WHERE (oexp2(1:NGI) > 300.0D0) !prevent floating point underflow as the second term becomes practically zero:
+ WHERE (oexp2(1:NGI) > 300.0D0) !prevent floating point underflow as the second term becomes practically zero:
Cnca(1:NGI,J) = cnac1(1:NGI,J)
ELSEWHERE
Cnca(1:NGI,J) = cnac1(1:NGI,J) +cnac2(1:NGI,J)*EXP(-oexp2(1:NGI))
@@ -240,7 +242,7 @@ SUBROUTINE LR_MR_activity()
CSca(1:NGI,J) = lrw1*omega2(1:NGI,J)*(Cnca(1:NGI,J)-cnac1(1:NGI,J))
ENDDO
- DO CONCURRENT (J = 1:NGI)
+ DO CONCURRENT (J = 1:NGI)
IF (SI2 > 250.0D0) THEN !prevent floating point underflow as the second term becomes practically zero:
Bkionc(1:NGN,J) = bnc(1:NGN,J)
Bkiona(1:NGN,J) = bna(1:NGN,J)
@@ -412,7 +414,7 @@ END SUBROUTINE LR_MR_activity
!****************************************************************************************
SUBROUTINE MRinteractcoeff()
- !Module variables:
+ !Module variables:
USE ModSystemProp, ONLY : NGN, NGI, NG, Ication, Ianion, Imaingroup, Ncation, isPEGsystem, &
& CatNr, AnNr, errorflagmix, frominpfile
@@ -589,7 +591,7 @@ SUBROUTINE MRinteractcoeff()
END SUBROUTINE MRinteractcoeff
!==========================================================================================================================
-
+
!****************************************************************************************
!* :: Purpose :: *
@@ -663,7 +665,17 @@ SUBROUTINE MRdata()
cn2TABAC(1,21) = -0.263257805913055D0
omega2TAB(1,21) = 1.31696686093901D0
TABhighestwtf(1,21) = 0.25D0
- TABKsp(1,21) = 3.25D0 ![molal basis] molal ion activity product at the solubility limit @ 298.15 K
+ TABKsp(1,21) = 3.25D0 ![molal basis] molal ion activity product at the solubility limit @ 298.15 K
+
+ !Li+ <--> HSO4- using an analogy approach with the (HSO4- <--> H+) parameters
+ bTABAC(1,8) = 2.15532299D-02
+ cTABAC(1,8) = 5.62965674D-01
+ omegaTAB(1,8) = 1.42442019D-01
+ cn1TABAC(1,8) = 0.0D0
+ cn2TABAC(1,8) = 7.03842440D-02
+ omega2TAB(1,8) = 7.14194282D-01
+ TABhighestwtf(1,8) = 0.82D0
+ TABKsp(1,8) = 1.0D28 ![molal basis] here just set to indicate suppressed solid formation [not true value]
!Na+ <-> Cl- 4P JUN 2007
bTABAC(2,2) = 5.374079546760321D-002
@@ -682,6 +694,15 @@ SUBROUTINE MRdata()
omega2TAB(2,3) = 2.20904992991060D0
TABhighestwtf(2,3) = 0.49D0
TABKsp(2,3) = 322.0D0 ![molal basis] molal ion activity product at the solubility limit @ 298 K
+
+ !Na+ <-> I-
+ bTABAC(2,4) = 0.1D0 !values to be determined...
+ cTABAC(2,4) = 0.1D0
+ cn1TABAC(2,4) = 0.0D0
+ cn2TABAC(2,4) = -0.1D0
+ omega2TAB(2,4) = 0.1D0
+ TABhighestwtf(2,4) = 0.6D0
+ TABKsp(2,4) = 999999.0D0 ![molal basis] , here just set to indicate suppressed solid formation, [not true value]!
!Na+ <-> NO3- 5P JUN 2007
bTABAC(2,5) = 1.164350105903698D-003
@@ -739,14 +760,23 @@ SUBROUTINE MRdata()
omegaTAB(3,5) = 0.357227451106132D0
TABhighestwtf(3,5) = 0.26D0
TABKsp(3,5) = 0.9225D0 ![molal basis] molal ion activity product at the solubility limit @ 298 K
-
- !K+ <-> SO4-- 4P JUN 2007 ; only few data points to constrain parameters!
- bTABAC(3,21) = 4.079155785591866D-003
- cTABAC(3,21) = -0.869936141934613D0
+
+ !K+ <-> HSO4- Revised May 2019
+ bTABAC(3,8) = 8.4433898901815352D-02
+ cTABAC(3,8) = -3.9301228443995739D-01
+ omegaTAB(3,8) = 0.8D0
+ cn1TABAC(3,8) = 0.0D0
+ cn2TABAC(3,8) = -9.0232965273929300D-01
+ omega2TAB(3,8) = 1.8692081240631140D+00
+ TABhighestwtf(3,8) = 0.90D0 ![not true value]
+
+ !K+ <-> SO4-- Revised May 2019
+ bTABAC(3,21) = 6.6714599718481538D-02
+ cTABAC(3,21) = -1.0844118398860187D+00
cn1TABAC(3,21) = 0.0D0
- cn2TABAC(3,21) = -9.224043798737591D-002
- omega2TAB(3,21) = 0.918742767314183D0
- TABhighestwtf(3,21) = 0.11D0
+ cn2TABAC(3,21) = -6.6511262694399914D-02
+ omega2TAB(3,21) = 5.8327209615736464D-01
+ TABhighestwtf(3,21) = 0.11D0
!NH4+ <-> Cl- 5P JUN 2007
bTABAC(4,2) = 1.520261933644218D-003
@@ -788,10 +818,10 @@ SUBROUTINE MRdata()
TABhighestwtf(4,21) = 0.78D0
TABKsp(4,21) = 1.03D0 ![molal basis]
- !NH4+ <-> HSO4- 5P FEB 2011
+ !NH4+ <--> HSO4- 5P FEB 2011
bTABAC(4,8) = 7.59734841D-03
cTABAC(4,8) = 1.43011653D-01
- omegaTAB(4,8) = 2.03954033D-01
+ omegaTAB(4,8) = 2.03954033D-01
cn1TABAC(4,8) = 0.0D0
cn2TABAC(4,8) = 6.31184288D-03
omega2TAB(4,8) = 8.25385786D-01
@@ -805,7 +835,7 @@ SUBROUTINE MRdata()
omega2TAB(5,2) = 0.504672383721131D0
TABhighestwtf(5,2) = 0.18D0
- !H+ <-> Br- 4P JUN 2007
+ !H+ <--> Br- 4P JUN 2007
bTABAC(5,3) = 0.120325353682262D0
cTABAC(5,3) = 0.444858839503779D0
cn1TABAC(5,3) = 0.0D0
@@ -813,7 +843,7 @@ SUBROUTINE MRdata()
omega2TAB(5,3) = 0.596775891883278D0
TABhighestwtf(5,3) = 0.20D0
- !H+ <-> NO3- 4P JUN 2007
+ !H+ <--> NO3- 4P JUN 2007
bTABAC(5,5) = 0.210637923183262D0
cTABAC(5,5) = 0.122694135055374D0
cn1TABAC(5,5) = 0.0D0
@@ -821,17 +851,17 @@ SUBROUTINE MRdata()
omega2TAB(5,5) = 1.67641985798924D0
TABhighestwtf(5,5) = 0.16D0
- !H2SO4 (-> HSO4- & H+) 5P FEB 2011
+ !H2SO4 (HSO4- <--> H+) 5P FEB 2011
bTABAC(5,8) = 2.15532299D-02
cTABAC(5,8) = 5.62965674D-01
- omegaTAB(5,8) = 1.42442019D-01 !!
+ omegaTAB(5,8) = 1.42442019D-01
cn1TABAC(5,8) = 0.0D0
cn2TABAC(5,8) = 7.03842440D-02
omega2TAB(5,8) = 7.14194282D-01
TABhighestwtf(5,8) = 0.82D0
TABKsp(5,8) = 1.0D28 ![molal basis] here just set to indicate suppressed solid formation of bisulfate and H+, [not true value]!
- !HSO4- (-> H+ & SO4--) 5P FEB 2011
+ !HSO4- ( H+ <--> SO4--) 5P FEB 2011
bTABAC(5,21) = 2.863428226D-01
cTABAC(5,21) = -5.996148056D0
omegaTAB(5,21) = 1.368612639D0
@@ -841,7 +871,7 @@ SUBROUTINE MRdata()
TABhighestwtf(5,21) = 0.82D0
TABKsp(5,21) = 1.0D28 ![molal basis] here just set to indicate suppressed solid formation of bisulfate, [not true value]!
- !Ca++ <-> NO3- 4P MAR 2006
+ !Ca++ <--> NO3- 4P MAR 2006
bTABAC(21,5) = 0.163281976992298D0
cTABAC(21,5) = 0.203681108454362D0
cn1TABAC(21,5) = 0.0D0
@@ -850,7 +880,7 @@ SUBROUTINE MRdata()
TABhighestwtf(21,5) = 0.50D0
TABKsp(21,5) = 2100.0D0 ![molal basis]
- !Ca++ <-> Cl- 4P JUN 2007
+ !Ca++ <--> Cl- 4P JUN 2007
bTABAC(21,2) = 0.104920450619807D0
cTABAC(21,2) = 0.866923035893255D0
cn1TABAC(21,2) = 0.0D0
@@ -858,13 +888,22 @@ SUBROUTINE MRdata()
omega2TAB(21,2) = 0.365747420374652D0
TABhighestwtf(21,2) = 0.53D0
- !Ca++ <-> Br- 4P JUL 2010
+ !Ca++ <--> Br- 4P JUL 2010
bTABAC(21,3) = 0.890929D0
cTABAC(21,3) = 6.101342D-002
cn1TABAC(21,3) = 0.0D0
cn2TABAC(21,3) = -0.238788D0
omega2TAB(21,3) = 0.762961D0
- TABhighestwtf(21,3) = 0.54D0
+ TABhighestwtf(21,3) = 0.54D0
+
+ !Ca++ <-> HSO4- AUG 2019
+ bTABAC(21,8) = 8.2098127298847412D-02
+ cTABAC(21,8) = -1.7166655056598561D+00
+ cn1TABAC(21,8) = 0.0D0
+ cn2TABAC(21,8) = 8.6726419889762707D-02
+ omega2TAB(21,8) = 4.6756611045939778D-01
+ TABhighestwtf(21,8) = 0.60D0 !Not True Value!
+ TABKsp(21,8) = 999999.0D0 ![molal basis]
! Ca++ <-> SO4-- 4P JAN 2011
bTABAC(21,21) = 1.2956723534082162D0
@@ -901,7 +940,16 @@ SUBROUTINE MRdata()
cn2TABAC(23,5) = -0.511836366410133D0
omega2TAB(23,5) = 1.44093984791712D0
TABhighestwtf(23,5) = 0.54D0
- TABKsp(23,5) = 30383.0D0 ![molal basis]
+ TABKsp(23,5) = 30383.0D0 ![molal basis]
+
+ !Mg++ <-> HSO4- JULY 2019
+ bTABAC(23,8) = 1.7141381771763914D-01
+ cTABAC(23,8) = 2.7138143502966559D+00
+ cn1TABAC(23,8) = 0.0D0
+ cn2TABAC(23,8) = 2.0404160088725892D-01
+ omega2TAB(23,8) = 5.0134675159022479D-01
+ TABhighestwtf(23,8) = 0.60D0 !Not True Value!
+ TABKsp(23,8) = 999999.0D0 ![molal basis]
! Mg++ <-> SO4-- 4P JUN 2007
bTABAC(23,21) = 0.122364180080768D0
diff --git a/FortranCode/ModSRunifac.f90 b/FortranCode/ModSRunifac.f90
index b8336a5..8e4d9ef 100644
--- a/FortranCode/ModSRunifac.f90
+++ b/FortranCode/ModSRunifac.f90
@@ -11,8 +11,8 @@
!* *
!* -> created: 2018 (based on non-module version from 2004) *
!* -> latest changes: 2018/05/25 *
-!* *
-!* :: License :: *
+!* *
+!* :: License :: *
!* This program is free software: you can redistribute it and/or modify it under the *
!* terms of the GNU General Public License as published by the Free Software *
!* Foundation, either version 3 of the License, or (at your option) any later *
@@ -36,7 +36,7 @@
!****************************************************************************************
MODULE ModSRunifac
-USE ModAIOMFACvar, ONLY : lastTK
+USE ModAIOMFACvar, ONLY : lastTK, eta0
USE ModSystemProp, ONLY : nd, nneutral, NG, NGN, Allsubs, ITABsr, Nmaingroups, topsubno, &
& maingrindexofsubgr, Imaingroup, isPEGsystem, calcviscosity, solvmixrefnd, COMPN
@@ -196,6 +196,7 @@ MODULE ModSRunifac
SUBROUTINE SRunifac(NK, T_K, X, XN, refreshgref, lnGaSR)
+ USE ModAIOMFACvar, ONLY : etamix, eta_cpn
IMPLICIT NONE
INTEGER(4),INTENT(IN) :: NK
@@ -206,24 +207,48 @@ SUBROUTINE SRunifac(NK, T_K, X, XN, refreshgref, lnGaSR)
!local variables:
INTEGER(4) :: J
REAL(8) :: RSS, expon
- REAL(8),DIMENSION(NK) :: lnGaC, lnGaR, phi
+ REAL(8),DIMENSION(NK) :: lnGaC, lnGaR, XieC, XieR, phi
!.......................................................
!Determine the reference values (lnGaRref, XieRref) for the residual part:
- IF (refreshgref .OR. solvmixrefnd) THEN
- CALL SRgref(NK, T_K, XN, refreshgref, lnGaRref)
+ IF (refreshgref .OR. solvmixrefnd .OR. calcviscosity) THEN
+ CALL SRgref(NK, T_K, XN, refreshgref, lnGaRref, XieRref)
ENDIF
!Residual part:
- CALL SRgres(X, lnGaR, NK)
+ CALL SRgres(X, lnGaR, NK, XieR)
!Combinatorial part:
- CALL SRgcomb(X, lnGaC, NK)
+ CALL SRgcomb(X, lnGaC, NK, XieC)
+ IF (calcviscosity) THEN !calculate the component's volume fraction phi
+ RSS = SUM(RS(:)*X(:))
+ phi(1:NK) = X(:)*RS(:)/RSS
+ ENDIF
+ etamix = 0.0D0
DO J = 1,NK
lnGaSR(J) = lnGaC(J) + lnGaR(J) - lnGaRref(J)
IF (J > nneutral) THEN !ion component, so normalize by reference state of infinite dilution in reference solvent
lnGaSR(J) = lnGaSR(J) - lnGaCinf(J)
ENDIF
+ IF (calcviscosity) THEN
+ !sum up liquid mixture viscosity contributions from different components:
+ eta_cpn(J) = XieC(J) +phi(J)*(XieR(J)-XieRref(J))
+ etamix = etamix + eta_cpn(J)
+ ENDIF
ENDDO
+ IF (calcviscosity) THEN
+ IF (etamix > -5.5D3) THEN !real calculation with valid parameters was performed
+ IF (etamix > 690.0D0) THEN !check and deal with floating point overflow issue
+ expon = DINT(etamix)
+ etamix = EXP(690.0D0 + 5.0D0*(etamix-expon))
+ ELSE
+ etamix = EXP(etamix) !"liquid" mixture viscosity in [Pa.s]
+ ENDIF
+ ELSE
+ etamix = -888.8D0 !indicate "not valid" calculation
+ ENDIF
+ ELSE
+ etamix = -999.9D0 !indicate "not calculated"
+ ENDIF
END SUBROUTINE SRunifac
!==========================================================================================================================
@@ -244,7 +269,10 @@ END SUBROUTINE SRunifac
!* -> latest changes: 2018/05/23 *
!* *
!****************************************************************************************
- SUBROUTINE SRgref(NK, T_K, XN, grefresh, lnGaRrefer)
+ SUBROUTINE SRgref(NK, T_K, XN, grefresh, lnGaRrefer, XieRrefer)
+
+ USE ModAIOMFACvar, ONLY : Tglass0, fragil ! public variable to hold pure component viscosity
+ USE Mod_PureViscosPar, ONLY : PureCompViscosity
IMPLICIT NONE
!interface variables:
@@ -252,20 +280,12 @@ SUBROUTINE SRgref(NK, T_K, XN, grefresh, lnGaRrefer)
REAL(8),INTENT(IN) :: T_K
REAL(8),DIMENSION(:),INTENT(IN) :: XN
LOGICAL(4),INTENT(IN) :: grefresh
- REAL(8),DIMENSION(:),INTENT(OUT) :: lnGaRrefer
- !interface to PureCompViscosity:
- INTERFACE
- PURE SUBROUTINE PureCompViscosity(compn, TempK, eta, iflag)
- INTEGER(4),INTENT(IN) :: compn
- REAL(8),INTENT(IN) :: TempK
- REAL(8),INTENT(OUT) :: eta
- INTEGER(4),INTENT(OUT) :: iflag
- END SUBROUTINE PureCompViscosity
- END INTERFACE
+ REAL(8),DIMENSION(:),INTENT(OUT) :: lnGaRrefer, XieRrefer
!local variables:
- INTEGER(4) :: I, iflag, cpn, K
+ INTEGER(4) :: I, iflag, K
REAL(8),PARAMETER :: T0 = 298.15D0 !room temperature as the reference temperature for the calculation of the PsiT array
- REAL(8),DIMENSION(NK) :: Xrefsp, lnGaRX
+ REAL(8),DIMENSION(NK) :: Xrefsp, lnGaRX, XieX
+ REAL(8) :: eta, Tglass, D ! pure component viscosity, glass transition temperature, fragility
!...................................................
!Check / set temperature-dependent parameters:
@@ -285,14 +305,16 @@ END SUBROUTINE PureCompViscosity
!reference residual contributions of pure components (and also hypothetical pure ions) are always 1.0 (i.e. for water)
Xrefsp = 0.0D0
lnGaRrefer = 0.0D0
+ XieRrefer = 0.0D0
DO I = 1,nneutral
!check if number of subgroups in component is greater than 1:
K = SUM(SRNY_dimflip(1:NGN,I))
IF (K > 1) THEN !more than one subgroup
Xrefsp(1:nneutral) = 0.0D0
Xrefsp(I) = 1.0D0 !set all other mole fractions but this (I) to 0.0. Thus reference state conditions of x(i) = 1
- CALL SRgres(Xrefsp, lnGaRX, NK)
+ CALL SRgres(Xrefsp, lnGaRX, NK, XieX)
lnGaRrefer(I) = lnGaRX(I) !this line is necessary since the other values (not only I) get overwritten in call to SRgres
+ XieRrefer(I) = XieX(I)
ENDIF
ENDDO
ENDIF
@@ -302,9 +324,47 @@ END SUBROUTINE PureCompViscosity
DO I = nneutral+1,NK
Xrefsp(1:nneutral) = XN(1:nneutral) !solvent
Xrefsp(nneutral+1:) = 0.0D0 !ions
- CALL SRgres(Xrefsp, lnGaRX, NK)
+ CALL SRgres(Xrefsp, lnGaRX, NK, XieX)
lnGaRrefer(I) = lnGaRX(I)
+ XieRrefer(I) = XieX(I)
+ ENDDO
+ ENDIF
+
+ !---- calculate pure component viscosity of non-electrolytes at given temperature
+ IF (calcviscosity) THEN
+ DO I = 1,nneutral
+ CALL PureCompViscosity(I, T_K, eta, iflag, Tglass, D) !calculate pure component dynamic viscosity eta
+ IF (iflag == 0) THEN
+ Tglass0(I) = Tglass
+ fragil(I) = D
+ eta0(I) = eta !eta in SI units of [Pa s]
+ ELSE
+ eta0(I) = -777.7D0 !to indicate a problem
+ SELECT CASE(nd)
+ CASE(500:799,2000:2500)
+ !do nothing for now....
+ CASE DEFAULT
+ !$OMP CRITICAL
+ WRITE(*,*) ""
+ WRITE(*,*) "WARNING: a pure compound viscosity for component",I,"could not be set due to missing parameters or a temperature outside the available bounds!"
+ WRITE(*,*) "Dataset is nd = ", nd
+ !READ(*,*)
+ !$OMP END CRITICAL
+ END SELECT
+ ENDIF
ENDDO
+ IF (NK > nneutral) THEN !calculate viscosity of (hypothetical) individual ions
+ !for electrolytes / ions:
+ !pure water as reference value
+ CALL PureCompViscosity(-1, T_K, eta, iflag, Tglass, D) !argument "-1" indicating water.
+ IF (iflag /= 0) THEN
+ eta = 1.0D-3 !eta in SI units of [Pa s]
+ ENDIF
+ !ions / electrolyte components:
+ eta0(nneutral+1:NK) = eta*100.0D0 !empirical value for dissolved ions (here set as 100 times the value of pure water viscosity)
+ ENDIF
+ ELSE
+ eta0(1:NK) = -999.9D0 !an impossible value to indicate uninitialized eta0
ENDIF
END SUBROUTINE SRgref
@@ -327,12 +387,12 @@ END SUBROUTINE SRgref
!* -> latest changes: 2018/05/18 *
!* *
!****************************************************************************************
- PURE SUBROUTINE SRgres(Xr, lnGaR, NK)
+ PURE SUBROUTINE SRgres(Xr, lnGaR, NK, XieR)
IMPLICIT NONE
!interface variables:
INTEGER(4),INTENT(IN) :: NK
- REAL(8),DIMENSION(:),INTENT(OUT) :: lnGaR
+ REAL(8),DIMENSION(:),INTENT(OUT) :: lnGaR, XieR
REAL(8),DIMENSION(:),INTENT(IN) :: Xr
!local variables:
INTEGER(4) :: I, k
@@ -359,6 +419,28 @@ PURE SUBROUTINE SRgres(Xr, lnGaR, NK)
lnGaR(I) = SUM(SRNY_dimflip(1:NG,I)*GAML(1:NG))
ENDDO
+ !--------------------------------------------------------
+ !Viscosity of mixture calculation; residual part.
+ !(electrolyte effects could be considered within the pure-component eta0 contribution from water as the main solvent of ions).
+ IF (calcviscosity) THEN
+ DO k = 1,NG !k
+ DO I = 1,NG !m
+ S4(I) = SUM(TH(1:NG)*PsiT(1:NG,I))
+ THmn(k,I) = TH(k)*PsiT_dimflip(I,k)/S4(I) !TH(k)*PsiT(k,I)/S4(I)
+ ENDDO
+ ENDDO
+ !sum up the number of subgroups and their contributions to XieR in compound I:
+ DO I = 1,NK
+ !calculate viscosity contributions by individual subgroups:
+ DO k = 1,NG
+ Xiek(k) = (-1.0D0)*(-Q(k)/R(k))*Nvis(k,I)*SUM(THmn(1:NG,k)*LOG(PsiT(1:NG,k)))
+ ENDDO
+ XieR(I) = SUM(SRNY_dimflip(1:NG,I)*Xiek(1:NG))
+ ENDDO
+ ELSE !a default value indicating "NO viscosity calc."
+ XieR = 0.0D0
+ ENDIF
+ !--------------------------------------------------------
END SUBROUTINE SRgres
!==========================================================================================================================
@@ -378,13 +460,13 @@ END SUBROUTINE SRgres
!* -> latest changes: 2018/05/18 *
!* *
!****************************************************************************************
- SUBROUTINE SRgcomb(X, lnGaC, NK)
+ SUBROUTINE SRgcomb(X, lnGaC, NK, XieC)
IMPLICIT NONE
!interface variables:
INTEGER(4),INTENT(IN) :: NK
REAL(8),DIMENSION(:),INTENT(IN) :: X
- REAL(8),DIMENSION(:),INTENT(OUT) :: lnGaC
+ REAL(8),DIMENSION(:),INTENT(OUT) :: lnGaC, XieC
!local variables:
INTEGER(4) :: nnp1
REAL(8) :: QSSref, RSSref, QSS, RSS, Xsolvtot, XLSref
@@ -417,6 +499,23 @@ SUBROUTINE SRgcomb(X, lnGaC, NK)
gcombrefresh = .false.
ENDIF
+ !--------------------------------------------------------
+ !Viscosity of mixture calculation; combinatorial part.
+ !(electrolyte effects could be considered within the pure-component eta0 contribution from water as the main solvent of ions).
+ IF (calcviscosity) THEN
+ phi = X*ViV !volume based component fraction
+ QiQ = QS(1:NK)/QSS
+ phiQ = X*QiQ
+ !undefined eta0 would contain negative values and cause error below
+ IF (ALL(eta0(1:NK) > 0.0D0)) THEN !pure component values available
+ XieC(1:NK) = (EXP(lnGaC(1:NK))*X(1:NK))*LOG(eta0(1:NK))
+ ELSE
+ XieC = -55555.5D0
+ ENDIF
+ ELSE !a default value indicating "NO viscosity calc."
+ XieC = -99999.9D0
+ ENDIF
+ !--------------------------------------------------------
END SUBROUTINE SRgcomb
!==========================================================================================================================
diff --git a/FortranCode/ModSubgroupProp.f90 b/FortranCode/ModSubgroupProp.f90
index 43e9023..b5c34a0 100644
--- a/FortranCode/ModSubgroupProp.f90
+++ b/FortranCode/ModSubgroupProp.f90
@@ -10,8 +10,8 @@
!* *
!* -> created: 2018 (based on non-module version from 2009) *
!* -> latest changes: 2018/05/22 *
-!* *
-!* :: License :: *
+!* *
+!* :: License :: *
!* This program is free software: you can redistribute it and/or modify it under the *
!* terms of the GNU General Public License as published by the Free Software *
!* Foundation, either version 3 of the License, or (at your option) any later *
@@ -29,6 +29,7 @@
!* - SUBROUTINE SubgroupNames *
!* - SUBROUTINE MaingroupNames *
!* - SUBROUTINE cpsubgrstring *
+!* - SUBROUTINE O2C_H2C_component *
!* - SUBROUTINE OtoCandHtoCmix *
!* *
!****************************************************************************************
@@ -540,52 +541,119 @@ END SUBROUTINE cpsubgrstring
!=================================================================================================================================
- !***********************************************************************************************************
- !* *
- !* Subroutine to calculate the average (organic) elemental O:C and H:C ratios of a given mixture. *
- !* *
- !* (c) Andi Zuend, Div. Chemistry & Chemical Engineering, California Institute of Technology, 2012 *
- !***********************************************************************************************************
+ !****************************************************************************************
+ !* :: Purpose :: *
+ !* Subroutine to calculate the sum of carbons, hydrogens, and oxygens in a *
+ !* component of a system containing organics. Also includes the calculation of O:C *
+ !* and H:C ratios of the pure components. *
+ !* *
+ !* :: Authors & Copyright :: *
+ !* Natalie Gervasi, Andi Zuend, *
+ !* Dept. Atmospheric and Oceanic Sciences, McGill University *
+ !* *
+ !* -> created: 2012 *
+ !* -> latest changes: 2018/09/14 *
+ !* *
+ !****************************************************************************************
+ PURE SUBROUTINE O2C_H2C_component(ind, compC, compH, compO, O2C, H2C)
+
+ USE ModSystemProp, ONLY : ITAB, SolvSubs, NGN
+
+ IMPLICIT NONE
+
+
+ !..................................
+ !interface input:
+ INTEGER(4),INTENT(IN) :: ind !the component index number in current mixture (e.g. index location of component in ITAB)
+ REAL(8),INTENT(OUT) :: compC, compH, compO, O2C, H2C
+ !local variables:
+ INTEGER(4) :: isub, k
+ !..................................
+
+ !(1) determine number of oxygen, hydrogen and carbon for a given compound
+ compO = 0.0D0
+ compH = 0.0D0
+ compC = 0.0D0
+ !loop over subgroups to count the O, H, and C atoms:
+
+ DO k = 1,NGN !loop over organic subgroups (SolvSubs excl. water)
+ isub = SolvSubs(k) !subgoup
+ IF (isub /= 16) THEN !exclude water (= subgroup 16) from the calculations
+
+ IF (ITAB(ind,isub) > 0) THEN !subgroup is present
+
+ compO = compO +REAL(ITAB(ind,isub)*subgO(isub), KIND=8)
+ compH = compH +REAL(ITAB(ind,isub)*subgH(isub), KIND=8)
+ compC = compC +REAL(ITAB(ind,isub)*subgC(isub), KIND=8)
+ ENDIF
+
+ ENDIF
+ ENDDO
+
+ !(2) compute O:C and H:C of this pure component:
+ IF (compC > 0.0D0) THEN
+ O2C = compO/compC
+ H2C = compH/compC
+ ELSE !O:C and H:C are undefined, labeled as negative numbers.
+ O2C = -77.77777D0
+ H2C = -77.77777D0
+ ENDIF
+
+ END SUBROUTINE O2C_H2C_component
+ !=================================================================================================================================
+
+
+ !****************************************************************************************
+ !* *
+ !* Subroutine to calculate the average (organic) elemental O:C and H:C ratios of a *
+ !* given mixture. *
+ !* *
+ !* :: Author & Copyright :: *
+ !* Andi Zuend, *
+ !* Dept. Chem. Engineering, California Institute of Technology (2009 - 2012) *
+ !* Dept. Atmospheric and Oceanic Sciences, McGill University *
+ !* *
+ !* -> created: 2012 *
+ !* -> latest changes: 2018/09/14 *
+ !* *
+ !****************************************************************************************
PURE SUBROUTINE OtoCandHtoCmix(n, x, OtoCorgmix, HtoCorgmix)
- USE ModSystemProp, ONLY : ITAB, nneutral, SolvSubs, NGN
+ USE ModSystemProp, ONLY : waterpresent
IMPLICIT NONE
!interface variables:
- INTEGER(4),INTENT(IN) :: n !number of species in mixture (potentially with electrolytes dissociated into ions)
- REAL(8),DIMENSION(n),INTENT(IN) :: x !mole fractions of species in mixture
+ INTEGER(4),INTENT(IN) :: n !number of species in mixture (potentially with electrolytes dissociated into ions)
+ REAL(8),DIMENSION(n),INTENT(IN) :: x !mole fractions of the species in the mixture
REAL(8),INTENT(OUT) :: OtoCorgmix, HtoCorgmix
- !..
!local variables:
- INTEGER(4) :: i, k, isub
- REAL(8) :: sumO, sumH, sumC, subamount
+ INTEGER(4) :: ind, istart
+ REAL(8) :: compO, compH, compC, sumO, sumH, sumC, O2C, H2C
!...........................................................
- !initialize:
- sumO = 0.0D0
- sumH = 0.0D0
+
sumC = 0.0D0
- !loop over compounds and subgroups to count the O, H, and C atoms
- !present in the mixture at given composition (mole fraction x):
- DO k = 1,NGN !loop over organic subgroups (SolvSubs excl. water)
- isub = SolvSubs(k) !subgoup
- IF (isub /= 16) THEN !exclude water (= subgroup 16) from the calculations
- DO i = 1,nneutral !loop over neutral components
- IF (ITAB(i,isub) > 0) THEN !subgroup is present
- subamount = REAL(ITAB(i,isub), KIND=8)*x(i)
- sumO = sumO +subamount*REAL(subgO(isub), KIND=8)
- sumH = sumH +subamount*REAL(subgH(isub), KIND=8)
- sumC = sumC +subamount*REAL(subgC(isub), KIND=8)
- ENDIF
- ENDDO
- ENDIF
+ sumH = 0.0D0
+ sumO = 0.0D0
+ !loop over organic components and sum up the contributions to total oxygen, carbon and hydrogen atoms for given mole fractions in mixture:
+ IF (waterpresent) THEN
+ istart = 2
+ ELSE
+ istart = 1
+ ENDIF
+ DO ind = istart,n
+ CALL O2C_H2C_component(ind, compC, compH, compO, O2C, H2C)
+ sumC = sumC + x(ind)*compC
+ sumH = sumH + x(ind)*compH
+ sumO = sumO + x(ind)*compO
ENDDO
- !calculate ratios from the sums:
+
+ !calculate organic O:C and H:C ratios for the present mixture:
IF (sumC > 0.0D0) THEN !the ratios are defined
OtoCorgmix = sumO/sumC
HtoCorgmix = sumH/sumC
ELSE
- OtoCorgmix = -77.77777D0 !indicate not defined O:C
- HtoCorgmix = -77.77777D0 !indicate not defined H:C
+ OtoCorgmix = -77.77777D0 !indicate O:C not defined
+ HtoCorgmix = -77.77777D0 !indicate H:C not defined
ENDIF
END SUBROUTINE OtoCandHtoCmix
diff --git a/FortranCode/ModSystemProp.f90 b/FortranCode/ModSystemProp.f90
index 7e433a4..fb32f6c 100644
--- a/FortranCode/ModSystemProp.f90
+++ b/FortranCode/ModSystemProp.f90
@@ -12,8 +12,8 @@
!* *
!* -> created: 2005 *
!* -> latest changes: 2018/05/26 *
-!* *
-!* :: License :: *
+!* *
+!* :: License :: *
!* This program is free software: you can redistribute it and/or modify it under the *
!* terms of the GNU General Public License as published by the Free Software *
!* Foundation, either version 3 of the License, or (at your option) any later *
@@ -69,7 +69,7 @@ MODULE ModSystemProp
REAL(8),DIMENSION(:),ALLOCATABLE :: K_el, nuestoich, SubGroupMW
REAL(8),DIMENSION(:),ALLOCATABLE :: Mmass !component molar mass in order of neutrals then electrolytes
!--
-CHARACTER(LEN=60),DIMENSION(:),ALLOCATABLE,PUBLIC :: cpname, compname, compnameTeX !component names in order of mixture components
+CHARACTER(LEN=200),DIMENSION(:),ALLOCATABLE,PUBLIC :: cpname, compname, compnameTeX !component names in order of mixture components
CHARACTER(LEN=16),DIMENSION(:),ALLOCATABLE,PUBLIC :: ionname, ionnameTeX
CHARACTER(LEN=3000),DIMENSION(:),ALLOCATABLE,PUBLIC :: compsubgroups, compsubgroupsTeX, compsubgroupsHTML
!--
@@ -84,7 +84,7 @@ MODULE SUBROUTINE SetSystem(ndi, datafromfile, ninp, cpnameinp, cpsubginp)
LOGICAL(4),INTENT(IN) :: datafromfile
INTEGER(4),INTENT(IN) :: ninp
!optional arguments at call:
- CHARACTER(LEN=60),DIMENSION(ninp),INTENT(IN), OPTIONAL :: cpnameinp
+ CHARACTER(LEN=200),DIMENSION(ninp),INTENT(IN), OPTIONAL :: cpnameinp
INTEGER(4),DIMENSION(ninp,topsubno),INTENT(IN), OPTIONAL :: cpsubginp
END SUBROUTINE SetSystem
!--
diff --git a/FortranCode/Mod_PureViscosPar.f90 b/FortranCode/Mod_PureViscosPar.f90
new file mode 100644
index 0000000..a23a90d
--- /dev/null
+++ b/FortranCode/Mod_PureViscosPar.f90
@@ -0,0 +1,199 @@
+!****************************************************************************************
+!* :: Purpose :: *
+!* Module providing access to the pure component viscosity parameter arrays. *
+!* *
+!* :: Author & Copyright :: *
+!* Andi Zuend, Natalie Gervasi *
+!* Dept. Atmospheric and Oceanic Sciences, McGill University *
+!* *
+!* -> created: 2012/09/07 *
+!* -> latest changes: 2019/08/23 *
+!* *
+!* :: List of subroutines and functions contained in this module: *
+!* -------------------------------------------------------------- *
+!* - SUBROUTINE DeRieux_Tno_Est *
+!* - SUBROUTINE VogelTemp *
+!* *
+!****************************************************************************************
+
+MODULE Mod_PureViscosPar
+
+IMPLICIT NONE
+!..................................
+!public parameter arrays:
+REAL(8),DIMENSION(1500,2),PUBLIC :: CorrelTrange !structure: (lower T limit, upper T limit) in [K]
+!..................................
+!=================================================================================================================================
+ CONTAINS
+!=================================================================================================================================
+
+ !****************************************************************************************
+ !* :: Purpose :: *
+ !* Calculation of Tg in [K] using the model from DeRieux et al. 2018 (ACP). *
+ !* [cP] = 0.01 [Poise] = 0.01 [g/(cm.s)] = 0.001 [Pa.s] = 0.001 [N.s/m^2] *
+ !* *
+ !* :: Authors & Copyright :: *
+ !* Natalie Gervasi, Andi Zuend, *
+ !* Dept. Atmospheric and Oceanic Sciences, McGill University *
+ !* *
+ !****************************************************************************************
+ PURE SUBROUTINE DeRieux_Tno_Est(ind, expTg, Tg)
+
+ USE ModSubgroupProp, ONLY : O2C_H2C_component
+
+ IMPLICIT NONE
+ !..................................
+ !interface variables:
+ INTEGER(4),INTENT(IN) :: ind !index no. of component in ITAB
+ REAL(8),INTENT(IN) :: expTg !the experimental Tg (if available)
+ REAL(8),INTENT(OUT) :: Tg !glass transition temp [K]
+ !local variables:
+ REAL(8) :: nCO, bC, bH, bCH, bO, bCO, sumC, sumH, sumO, OtoC, HtoC
+ !..................................
+
+ CALL O2C_H2C_component(ind, sumC, sumH, sumO, OtoC, HtoC)
+
+ ! Tg is either calculated or an experimentally determined value is used
+ IF (expTg > -999.0D0) THEN ! use experimental value
+ Tg = expTg
+ ELSE
+ !choose the constants for the Tg model based on if the compound has oxygens or not
+ IF (sumO < 1.0D0) THEN
+ nCO = 1.96D0
+ bC = 61.99D0
+ bH = -113.33D0
+ bCH = 28.74D0
+ Tg = bC*(nCO + LOG(sumC)) + bH*LOG(sumH) + bCH*LOG(sumC)*LOG(sumH)
+ ELSE
+ nCO = 12.13D0
+ bC = 10.95D0
+ bH = -41.82D0
+ bCH = 21.61D0
+ bO = 118.96D0
+ bCO = -24.38D0
+ ! Shiraiwa et al. group model used to calculate Tg (DeRieux et al. 2018), eqn (2)
+ Tg = bC*(nCO + LOG(sumC)) + bH*LOG(sumH) + bCH*LOG(sumC)*LOG(sumH) + bO*LOG(sumO) + bCO*LOG(sumC)*LOG(sumO)
+ ENDIF
+ ENDIF
+
+ END SUBROUTINE DeRieux_Tno_Est
+ !=================================================================================================================================
+
+
+ !****************************************************************************************
+ !* :: Purpose :: *
+ !* Calculation of Tno in [K] using the Vogel-Tamman-Fulcher equation in the form *
+ !* used by DeRieux et al. 2018. Fragility (D) is determined by the temperature of the *
+ !* run. *
+ !* *
+ !* :: Authors & Copyright :: *
+ !* Natalie Gervasi, Andi Zuend, *
+ !* Dept. Atmospheric and Oceanic Sciences, McGill University *
+ !* *
+ !****************************************************************************************
+ PURE SUBROUTINE VogelTemp(Tg, TempK, D, Tno)
+
+ IMPLICIT NONE
+
+ !interface
+ REAL(8), INTENT(IN) :: Tg, TempK ! glass transition temperature, run temperature [K]
+ REAL(8), INTENT(OUT) :: Tno, D ! vogel temperature, fragility
+ !........................
+
+ IF (TempK < Tg) THEN ! pure component substance below Tg that was once fragile will behave more strongly
+ D = 30.0D0
+ ELSE
+ D = 10.0D0 ! reasonable guess for organics (Shiraiwa et al., 2017)
+ ENDIF
+ Tno = (39.17D0*Tg)/(D + 39.17D0) ! (DeRieux et al. 2018) eqn (7)
+
+ END SUBROUTINE VogelTemp
+ !=================================================================================================================================
+
+
+ !****************************************************************************************
+ !* :: Purpose :: *
+ !* Calculation of pure component dynamic viscosity (eta) in Pascal seconds [Pa.s] at *
+ !* a given temperature T [K]. *
+ !* [cP] = 0.01 [Poise] = 0.01 [g/(cm.s)] = 0.001 [Pa.s] = 0.001 [N.s/m^2] *
+ !* *
+ !* :: Author & Copyright :: *
+ !* Andi Zuend, Natalie Gervasi *
+ !* Dept. Atmospheric and Oceanic Sciences, McGill University *
+ !* *
+ !* -> created: 2012/09/07 *
+ !* -> latest changes: 2019/03/16 *
+ !* *
+ !****************************************************************************************
+
+ SUBROUTINE PureCompViscosity(ind, TempK, eta, iflag, Tglass, fragility)
+
+ USE ModSystemProp, ONLY : CompN, ITAB
+
+ IMPLICIT NONE
+ !..................................
+ !interface variables:
+ INTEGER(4),INTENT(IN) :: ind
+ REAL(8),INTENT(IN) :: TempK
+ REAL(8),INTENT(OUT) :: Tglass, eta, fragility
+ INTEGER(4),INTENT(OUT) :: iflag
+ !local variables:
+ INTEGER(4) :: cpn, equationNo
+ REAL(8) :: a, b, c, d, e, Tg, Tvog
+ !..................................
+
+ !initialize:
+ IF (ind > 0) THEN
+ IF (ITAB(ind,16) > 0) THEN !is water
+ cpn = 401
+ ELSE
+ cpn = compN(ind) !organic or water component ID (when defined)
+ ENDIF
+ ELSE IF (ind == -1) THEN
+ cpn = 401 !water
+ ENDIF
+
+ SELECT CASE(cpn)
+ CASE(401) !water
+ equationNo = 1
+ a = 225.66D0
+ b = 1.3788D-4
+ c = -1.6433D0
+ d = 0.0D0
+ e = 0.0D0
+ CorrelTrange(cpn,1) = 230.0D0
+ CorrelTrange(cpn,2) = 495.0D0
+ CASE DEFAULT !for system input of organics from file
+ equationNo = 2
+ a = -999.0D0
+ b = -999.0D0
+ c = -999.0D0
+ d = -999.0D0
+ e = -999.0D0
+ CorrelTrange(cpn,1) = 1.0D0
+ CorrelTrange(cpn,2) = 1000.0D0
+ END SELECT
+
+ !compare requested temperature with temperature range of parameterization:
+ IF (TempK >= CorrelTrange(cpn,1) .AND. TempK <= CorrelTrange(cpn,2)) THEN !valid
+ iflag = 0 !0 means no errors
+ ! decide apppropriate value in Tg range
+ CALL DeRieux_Tno_Est(ind, a, Tg)
+ Tglass = Tg
+ CALL VogelTemp(Tglass, TempK, fragility, Tvog)
+ !compute eta with given equation for the component
+ SELECT CASE(equationNo)
+ CASE(1)
+ eta = b*((TempK/a) - 1)**c
+ CASE(2) !Vogel-Fulcher-Tammann (VFT) using DeRieux et al. (2018) constants and DeRieux Tg
+ eta = 10.0D0**(-5.0D0 + 0.434D0*(fragility*Tvog/(TempK - Tvog))) !DeRieux et al. 2018, eqn (6)
+ END SELECT
+ iflag = 0 !0 means no errors
+ ELSE
+ iflag = 1 !1 = outside valid temperature range!
+ ENDIF
+
+ END SUBROUTINE PureCompViscosity
+ !=================================================================================================================================
+
+END MODULE Mod_PureViscosPar
\ No newline at end of file
diff --git a/FortranCode/OutputHTML.f90 b/FortranCode/OutputHTML.f90
new file mode 100644
index 0000000..410faf4
--- /dev/null
+++ b/FortranCode/OutputHTML.f90
@@ -0,0 +1,281 @@
+!****************************************************************************************
+!* :: Purpose :: *
+!* Create an output file in HTML markup format with an overall mixture data header *
+!* and individual tables for all components / species (in case of ions). *
+!* This output format is used for showning the results on the AIOMFAC-web webpage. *
+!* *
+!* :: Author & Copyright :: *
+!* Andi Zuend, (andi.zuend@gmail.com) *
+!* Dept. Atmospheric and Oceanic Sciences, McGill University (2013 - present) *
+!* *
+!* -> created: 2011 *
+!* -> latest changes: 2018/08/10 *
+!* *
+!* :: License :: *
+!* This program is free software: you can redistribute it and/or modify it under the *
+!* terms of the GNU General Public License as published by the Free Software *
+!* Foundation, either version 3 of the License, or (at your option) any later *
+!* version. *
+!* The AIOMFAC model code 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 General Public License for more *
+!* details. *
+!* You should have received a copy of the GNU General Public License along with this *
+!* program. If not, see . *
+!* *
+!****************************************************************************************
+SUBROUTINE OutputHTML(fname, VersionNo, cpnameinp, nspecmax, npoints, watercompno, T_K, px, out_data, out_viscdata)
+
+!module variables:
+USE ModSystemProp, ONLY : compname, compsubgroupsHTML, NGS, NKNpNGS, ninput, nneutral
+USE ModSubgroupProp, ONLY : subgrnameHTML
+
+IMPLICIT NONE
+!interface variables:
+CHARACTER(LEN=*),INTENT(IN) :: fname
+CHARACTER(LEN=*),INTENT(IN) :: VersionNo
+CHARACTER(LEN=*),DIMENSION(nspecmax),INTENT(IN) :: cpnameinp !list of assigned component names (from input file)
+INTEGER(4),INTENT(IN) :: nspecmax, npoints, watercompno
+REAL(8),DIMENSION(npoints),INTENT(IN) :: T_K
+INTEGER(4),DIMENSION(nspecmax),INTENT(INOUT) :: px
+REAL(8),DIMENSION(7,npoints,NKNpNGS),INTENT(IN) :: out_data
+REAL(8),DIMENSION(3,npoints),INTENT(IN) :: out_viscdata
+!--
+!local variables:
+CHARACTER(LEN=:),ALLOCATABLE :: cn
+CHARACTER(LEN=5) :: tlen
+CHARACTER(LEN=50) :: subntxt
+CHARACTER(LEN=150) :: cnformat, tformat, txtn
+CHARACTER(LEN=400) :: outtxtleft
+CHARACTER(LEN=3000) :: txtsubs, mixturestring
+INTEGER(4) :: i, k, pointi, qty, unitx
+REAL(8) :: RH
+!...................................................................................
+
+OPEN (NEWUNIT = unitx, FILE = fname, STATUS = "UNKNOWN")
+outtxtleft = ADJUSTL("AIOMFAC-web, version "//VersionNo//"
")
+WRITE(unitx,'(A400)') outtxtleft
+!create a character string of the mixture as a series of its components and add links to the components:
+mixturestring = ''//TRIM(ADJUSTL(compname(1)))//'' !first component
+!loop over all further components / species:
+DO i = 2,nspecmax
+ IF (INT(out_data(6,px(i),i)) == 0) THEN !neutral component
+ txtn = ''//TRIM(ADJUSTL(compname(i)))//''
+ ELSE !ion (with its own link)
+ txtn = ''//TRIM(ADJUSTL(subgrnameHTML(INT(out_data(6,px(i),i) ))))//''
+ ENDIF
+ mixturestring = TRIM(mixturestring)//' + '//TRIM(ADJUSTL(txtn))
+ENDDO
+mixturestring = ADJUSTL(TRIM(mixturestring))
+WRITE(tlen,'(I5.5)') LEN_TRIM(mixturestring)
+tformat = '(A24, A'//tlen//')' !dynamic format specifier
+WRITE(unitx, tformat) "Mixture name: ", ADJUSTL(mixturestring)
+WRITE(unitx,'(A45, I2.2)') "
Number of independent input components: ", ninput
+WRITE(unitx,'(A45, I2.2)') "
Number of different neutral components: ", nneutral
+WRITE(unitx,'(A45, I2.2)') "
Number of different inorganic ions : ", NGS
+WRITE(unitx,*) "
"
+outtxtleft = ADJUSTL(' The AIOMFAC output is tabulated for each component/species individually below.
')
+WRITE(unitx,'(A400)') outtxtleft
+WRITE(unitx,*) '
'
+!--
+outtxtleft = ADJUSTL('')
+WRITE(unitx,'(A400)') outtxtleft
+outtxtleft = ADJUSTL(' Table key ')
+WRITE(unitx,'(A400)') outtxtleft
+outtxtleft = ADJUSTL('| no. | : | composition point number; |
')
+WRITE(unitx,'(A400)') outtxtleft
+outtxtleft = ADJUSTL('| T [K] | : | absolute temperature; |
')
+WRITE(unitx,'(A400)') outtxtleft
+outtxtleft = ADJUSTL('| RH [%] | : | relative humidity in equilibrium with the liquid mixture (bulk conditions); |
')
+WRITE(unitx,'(A400)') outtxtleft
+outtxtleft = ADJUSTL('| w(j) [-] | : | weight fraction (mass fraction) of species "j"; |
')
+WRITE(unitx,'(A400)') outtxtleft
+outtxtleft = ADJUSTL('| x_i(j) [-] | : | mole fraction of species "j", calculated on the basis of completely dissociated inorganic ions; exception: the partial dissociation of bisulfate (HSO4- ↔ H+ + SO42-) is explicitly considered when present in the mixture; |
')
+WRITE(unitx,'(A400)') outtxtleft
+outtxtleft = ADJUSTL('| m_i(j) [mol/kg] | : | molality of species "j" [mol(j)/(kg solvent mixture)], where "solvent mixture" refers to the electrolyte-free mixture (water + organics); |
')
+WRITE(unitx,'(A400)') outtxtleft
+outtxtleft = ADJUSTL('| a_coeff_x(j) [-] | : | activity coefficient of "j", defined on mole fraction basis (used for non-ionic components) with pure (liquid) component "j" reference state; |
')
+WRITE(unitx,'(A400)') outtxtleft
+outtxtleft = ADJUSTL('| a_coeff_m(j) [-] | : | activity coefficient of "j", defined on molality basis (used for inorg. ions) with reference state of infinite dilution of "j" in pure water; |
')
+WRITE(unitx,'(A400)') outtxtleft
+outtxtleft = ADJUSTL('| a_x(j) [-] | : | activity of "j", defined on mole fraction basis (pure component reference state); |
')
+WRITE(unitx,'(A400)') outtxtleft
+outtxtleft = ADJUSTL('| a_m(j) [-] | : | activity of "j", defined on molality basis (used for inorg. ions) with reference state of infinite dilution of "j" in pure water; |
')
+WRITE(unitx,'(A400)') outtxtleft
+outtxtleft = ADJUSTL('| log10(η/[Pa.s]) | : | base-10 log of the dynamic viscosity of the mixture; |
')
+WRITE(unitx,'(A400)') outtxtleft
+outtxtleft = ADJUSTL('| log10(η sens./[Pa.s]) | : | base-10 log of the estimated sensitivity of the dynamic viscosity of the mixture; see details here; |
')
+WRITE(unitx,'(A400)') outtxtleft
+outtxtleft = ADJUSTL('| flag | : | error/warning flag, a non-zero value (error/warning number) indicates that a numerical issue or a warning occurred at this data point, suggesting evaluation with caution (warnings) or exclusion (errors) of this data point; see also the List_of_Error_and_Warning_Codes.pdf. |
')
+WRITE(unitx,'(A400)') outtxtleft
+outtxtleft = ADJUSTL('
')
+WRITE(unitx,'(A400)') outtxtleft
+WRITE(unitx,*) '
'
+WRITE(unitx,*) '
'
+
+!--
+!data table for mixture viscosity output
+WRITE(unitx,*) '
'
+!write component properties table above data table:
+outtxtleft = ADJUSTL('')
+WRITE(unitx,'(A400)') outtxtleft
+WRITE(unitx,*) ''
+WRITE(unitx,'(A61)') "| Properties of this phase: mixture viscosity |
"
+WRITE(unitx,*) ''
+WRITE(unitx,*) '
'
+!write table column headers:
+outtxtleft = ADJUSTL('')
+WRITE(unitx,'(A400)') outtxtleft
+WRITE(unitx,*) ''
+txtsubs = ADJUSTL("| no. | | T [K] | | RH [%] | | log10(η/[Pa.s]) | | log10(η sens./[Pa.s]) | | flag |
|---|
")
+WRITE(tlen,'(I0)') LEN_TRIM(txtsubs)
+tformat = '(A'//tlen//')'
+WRITE(unitx,tformat) ADJUSTL(txtsubs)
+WRITE(unitx,*) ''
+WRITE(unitx,*) ''
+!--
+!write data to table:
+DO pointi = 1,npoints !loop over composition points
+ IF (watercompno > 0) THEN !water is present in mixture
+ RH = out_data(5,pointi,watercompno)*100.0D0 !RH in %
+ IF (RH > 1000.0D0 .OR. RH < 0.0D0) THEN
+ RH = -99.99D0
+ ENDIF
+ ELSE
+ RH = 0.0D0
+ ENDIF
+ tformat = '(A8,I5.3," | ",F7.2," | | ",F7.2," | | ",2(ES12.5," | | "),I2,A10)'
+ WRITE(unitx, tformat) " | | ", pointi, T_K(pointi), RH, out_viscdata(1,pointi), out_viscdata(2,pointi), INT(out_viscdata(3,pointi)), " |
"
+ENDDO !pointi
+WRITE(unitx,*) ''
+WRITE(unitx,*) '
'
+WRITE(unitx,*) '↑ Top'
+WRITE(unitx,*) "
"
+WRITE(unitx,*) "
"
+!--
+k = MAX(2, CEILING(LOG10(REAL(nspecmax))) )
+WRITE(tlen,'(I0)') k
+cnformat = "(I"//TRIM(tlen)//"."//TRIM(tlen)//")" !variable format specifier
+ALLOCATE(CHARACTER(LEN=k) :: cn)
+
+!write individual data tables for each component / ionic species:
+DO i = 1,nspecmax
+ WRITE(cn,cnformat) i !component / species number as character string
+ IF (INT(out_data(6,px(i),i)) == 0) THEN !neutral component
+ WRITE(unitx,*) '
'
+ txtn = ADJUSTL('')
+ WRITE(unitx,'(A150)') txtn
+ !write component properties table above data table:
+ outtxtleft = ADJUSTL('')
+ WRITE(unitx,'(A400)') outtxtleft
+ WRITE(unitx,*) ''
+ WRITE(unitx,'(A73,I2.2,A10)') "| Mixture's component, # | : ", i ," |
"
+ txtn = TRIM(ADJUSTL(compname(i)))
+ txtn = TRIM(txtn)//""
+ WRITE(tlen,'(I5.5)') LEN_TRIM(txtn)
+ tformat = '(A66, A'//tlen//')'
+ WRITE(unitx, tformat) "| Component's name | : ", ADJUSTL(txtn)
+ txtsubs = ''//TRIM(compsubgroupsHTML(i))//" |
"
+ WRITE(tlen,'(I5.5)') LEN_TRIM(txtsubs)
+ tformat = '(A71, A'//tlen//')'
+ WRITE(unitx, tformat) "| Component's subgroups | : ", ADJUSTL(txtsubs)
+ WRITE(unitx,*) ' |
'
+ WRITE(unitx,*) '
'
+ !write table column headers:
+ outtxtleft = ADJUSTL('')
+ WRITE(unitx,'(A400)') outtxtleft
+ WRITE(unitx,*) ''
+ outtxtleft = ADJUSTL("| no. | | T [K] | | RH [%] | | w("//cn//") | | x_i("//cn//") | | m_i("//cn//") | | a_coeff_x("//cn//") | | a_x("//cn//") | | flag |
|---|
")
+ WRITE(unitx,'(A400)') outtxtleft
+ !--
+ ELSE IF (INT(out_data(6,px(i),i)) < 240) THEN !cation
+ WRITE(unitx,*) '
'
+ txtn = ''
+ WRITE(unitx,'(A150)') txtn !link target
+ !---
+ !write component properties table above data table:
+ outtxtleft = ADJUSTL('')
+ WRITE(unitx,'(A400)') outtxtleft
+ WRITE(unitx,*) ''
+ WRITE(unitx,'(A71,I2.2,A10)') "| Mixture's species, # | : ", i ," |
"
+ subntxt = TRIM(ADJUSTL(subgrnameHTML(INT(out_data(6,px(i),i)))))
+ qty = LEN_TRIM(subntxt)
+ txtn = ADJUSTL(subntxt(2:qty-1)) !to print the ion name without the enclosing parathesis ()
+ txtn = ''//TRIM(txtn)//""
+ WRITE(tlen,'(I5.5)') LEN_TRIM(txtn)
+ tformat = '(A63, A'//tlen//')'
+ WRITE(unitx, tformat) "| Cation's name | : ", ADJUSTL(txtn)
+ txtsubs = ''//TRIM(subntxt)//" |
"
+ WRITE(tlen,'(I5.5)') LEN_TRIM(txtsubs)
+ tformat = '(A68, A'//tlen//')'
+ WRITE(unitx, tformat) "| Cation's subgroups | : ", ADJUSTL(txtsubs)
+ WRITE(unitx,*) ' |
'
+ WRITE(unitx,*) '
'
+ !write data table column headers:
+ outtxtleft = ''
+ WRITE(tlen,'(I5.5)') LEN_TRIM(outtxtleft)
+ tformat = '(A'//tlen//')'
+ WRITE(unitx, tformat) ADJUSTL(outtxtleft)
+ WRITE(unitx,*) ''
+ outtxtleft = ADJUSTL('| no. | | T [K] | | RH [%] | | w('//cn//') | | x_i('//cn//') | | m_i('//cn//') | | a_coeff_m('//cn//') | | a_m('//cn//') | | flag |
|---|
')
+ WRITE(unitx,'(A400)') outtxtleft
+ !--
+ ELSE IF (INT(out_data(6,px(i),i)) > 240) THEN !anion
+ WRITE(unitx,*) '
'
+ txtn = ''
+ WRITE(unitx,'(A150)') txtn
+ !---
+ !write component properties table above data table:
+ outtxtleft = ADJUSTL('')
+ WRITE(unitx,'(A400)') outtxtleft
+ WRITE(unitx,*) ''
+ WRITE(unitx,'(A71,I2.2,A10)') "| Mixture's species, # | : ", i ," |
"
+ subntxt = TRIM(ADJUSTL(subgrnameHTML(INT(out_data(6,px(i),i)))))
+ qty = LEN_TRIM(subntxt)
+ txtn = ADJUSTL(subntxt(2:qty-1)) !to print the ion name without the enclosing parathesis ()
+ txtn = ''//TRIM(txtn)//""
+ WRITE(tlen,'(I5.5)') LEN_TRIM(txtn)
+ tformat = '(A63, A'//tlen//')'
+ WRITE(unitx, tformat) "| Anion's name | : ", ADJUSTL(txtn)
+ txtsubs = ''//TRIM(subntxt)//" |
"
+ WRITE(tlen,'(I5.5)') LEN_TRIM(txtsubs)
+ tformat = '(A68, A'//tlen//')'
+ WRITE(unitx, tformat) "| Anion's subgroups | : ", ADJUSTL(txtsubs)
+ WRITE(unitx,*) ' |
'
+ WRITE(unitx,*) '
'
+ !write data table column headers:
+ outtxtleft = ADJUSTL('')
+ WRITE(unitx,'(A400)') outtxtleft
+ WRITE(unitx,*) ''
+ outtxtleft = ADJUSTL("| no. | | T [K] | | RH [%] | | w("//cn//") | | x_i("//cn//") | | m_i("//cn//") | | a_coeff_m("//cn//") | | a_m("//cn//") | | flag |
|---|
")
+ WRITE(unitx,'(A400)') outtxtleft
+ !--
+ ELSE
+ !error
+ ENDIF
+ !--
+ WRITE(unitx,*) ''
+ WRITE(unitx,*) ''
+ !write data to table:
+ DO pointi = 1,npoints !loop over composition points
+ IF (watercompno > 0) THEN !water is present in mixture
+ RH = out_data(5,pointi,watercompno)*100.0D0 !RH in %
+ IF (RH > 1000.0D0 .OR. RH < 0.0D0) THEN
+ RH = -99.99D0
+ ENDIF
+ ELSE
+ RH = 0.0D0
+ ENDIF
+ tformat = '(A8,I5.3," | ",F7.2," | | ",F7.2," | | ",5(ES12.5," | | "),I2,A10)'
+ WRITE(unitx, tformat) " | | ", pointi, T_K(pointi), RH, (out_data(k,pointi,i), k = 1,5), INT(out_data(7,pointi,i)), " |
"
+ ENDDO !pointi
+ WRITE(unitx,*) ''
+ WRITE(unitx,*) '
'
+ WRITE(unitx,*) '↑ Top'
+ WRITE(unitx,*) "
"
+ WRITE(unitx,*) "
"
+ENDDO
+
+CLOSE(unitx)
+
+END SUBROUTINE OutputHTML
\ No newline at end of file
diff --git a/FortranCode/OutputTXT.f90 b/FortranCode/OutputTXT.f90
index 170f104..7b3cac8 100644
--- a/FortranCode/OutputTXT.f90
+++ b/FortranCode/OutputTXT.f90
@@ -23,7 +23,7 @@
!* program. If not, see . *
!* *
!****************************************************************************************
-SUBROUTINE OutputTXT(fname, VersionNo, cpnameinp, nspecmax, npoints, watercompno, T_K, px, out_data)
+SUBROUTINE OutputTXT(fname, VersionNo, cpnameinp, nspecmax, npoints, watercompno, T_K, px, out_data, out_viscdata)
!module variables:
USE ModSystemProp, ONLY : compname, compsubgroups, compsubgroupsTeX, NGS, NKNpNGS, ninput, nneutral
@@ -33,17 +33,18 @@ SUBROUTINE OutputTXT(fname, VersionNo, cpnameinp, nspecmax, npoints, watercompno
!interface variables:
CHARACTER(LEN=*),INTENT(IN) :: fname
CHARACTER(LEN=*),INTENT(IN) :: VersionNo
-CHARACTER(LEN=60),DIMENSION(nspecmax),INTENT(IN) :: cpnameinp !list of assigned component names (from input file)
+CHARACTER(LEN=*),DIMENSION(nspecmax),INTENT(IN) :: cpnameinp !list of assigned component names (from input file)
INTEGER(4),INTENT(IN) :: nspecmax, npoints, watercompno
REAL(8),DIMENSION(npoints),INTENT(IN) :: T_K
INTEGER(4),DIMENSION(nspecmax),INTENT(INOUT) :: px
REAL(8),DIMENSION(7,npoints,NKNpNGS),INTENT(IN) :: out_data
+REAL(8),DIMENSION(3,npoints),INTENT(IN) :: out_viscdata
!--
!local variables:
-CHARACTER(LEN=2) :: cn
+CHARACTER(LEN=:),ALLOCATABLE :: cn
CHARACTER(LEN=5) :: tlen
CHARACTER(LEN=50) :: subntxt
-CHARACTER(LEN=150) :: horizline, tablehead, tformat, txtn
+CHARACTER(LEN=150) :: cnformat, horizline, tablehead, tformat, txtn
CHARACTER(LEN=3000) :: txtsubs, mixturestring
INTEGER(4) :: i, k, pointi, qty, unitx
REAL(8) :: RH
@@ -100,6 +101,9 @@ SUBROUTINE OutputTXT(fname, VersionNo, cpnameinp, nspecmax, npoints, watercompno
WRITE(unitx,'(A100)') 'a_x(j) [-] : activity of "j", defined on mole fraction basis (pure component reference state);'
WRITE(unitx,'(A100)') 'a_m(j) [-] : activity of "j", defined on molality basis (used for inorg. ions) with reference '
WRITE(unitx,'(A100)') ' state of infinite dilution of "j" in pure water; '
+WRITE(unitx,'(A100)') 'log_10(eta/[Pa.s]): base-10 log of the dynamic viscosity of the mixture; '
+WRITE(unitx,'(A100)') 'log_10(s_eta/[Pa.s]): base-10 log of the estimated sensitivity of the dynamic viscosity; see details'
+WRITE(unitx,'(A100)') ' on the "Hints & Examples" webpage; '
WRITE(unitx,'(A100)') 'flag : error/warning flag, a non-zero value (error/warning number) indicates that a '
WRITE(unitx,'(A100)') ' numerical issue or a warning occurred at this data point, suggesting evaluation '
WRITE(unitx,'(A100)') ' with caution (warnings) or exclusion (errors) of this data point. '
@@ -108,9 +112,36 @@ SUBROUTINE OutputTXT(fname, VersionNo, cpnameinp, nspecmax, npoints, watercompno
WRITE(unitx,*) ''
!--
horizline = "-----------------------------------------------------------------------------------------------------------"
+!--
+!data table for mixture viscosity output
+WRITE(unitx,'(A44)') "Properties of this phase: mixture viscosity"
+!write table column headers:
+WRITE(unitx,'(A107)') ADJUSTL(horizline)
+tablehead = "no. T_[K] RH_[%] log_10(eta/[Pa.s]) log_10(s_eta/[Pa.s]) flag "
+WRITE(unitx,'(2X, A105)') ADJUSTL(tablehead)
+!write data to viscosity table
+DO pointi = 1,npoints !loop over composition points
+ IF (watercompno > 0) THEN
+ RH = out_data(5,pointi,watercompno)*100.0D0 !RH in %
+ IF (RH > 1000.0D0 .OR. RH < 0.0D0) THEN
+ RH = -99.99D0
+ ENDIF
+ ELSE
+ RH = 0.0D0
+ ENDIF
+ WRITE(unitx,'(I5.3,2X,F7.2,2X,F7.2,3X,2(ES12.5,10X),3X,I2)') pointi, T_K(pointi), RH, out_viscdata(1,pointi), out_viscdata(2,pointi), INT(out_viscdata(3,pointi))
+ENDDO !pointi
+WRITE(unitx,'(A107)') ADJUSTL(horizline)
+WRITE(unitx,*) ""
+!--
+k = MAX(2, CEILING(LOG10(REAL(nspecmax))) )
+WRITE(tlen,'(I0)') k
+cnformat = "(I"//TRIM(tlen)//"."//TRIM(tlen)//")" !variable format specifier
+ALLOCATE(CHARACTER(LEN=k) :: cn)
+
!write individual data tables for each component / ionic species:
DO i = 1,nspecmax
- WRITE(cn,'(I2.2)') i !component / species number as character string
+ WRITE(cn,cnformat) i !component / species number as character string
WRITE(unitx,*) ""
!distinguish between neutral components and ionic species:
IF (INT(out_data(6,px(i),i)) == 0) THEN !neutral component
@@ -199,6 +230,8 @@ SUBROUTINE OutputTXT(fname, VersionNo, cpnameinp, nspecmax, npoints, watercompno
WRITE(unitx,*) ""
ENDDO
WRITE(unitx,*) ""
+!--
+WRITE(unitx,*) ""
WRITE(unitx,'(A107)') "==========================================================================================================="
CLOSE(unitx)
diff --git a/FortranCode/ReadInputFile.f90 b/FortranCode/ReadInputFile.f90
new file mode 100644
index 0000000..f4eaa32
--- /dev/null
+++ b/FortranCode/ReadInputFile.f90
@@ -0,0 +1,316 @@
+!****************************************************************************************
+!* :: Purpose :: *
+!* Read an input text file with information about the subgroups of each system *
+!* component and the compositions and temperatures for specific mixture calculations. *
+!* *
+!* :: Author & Copyright :: *
+!* Andi Zuend, (andi.zuend@gmail.com) *
+!* Dept. Atmospheric and Oceanic Sciences, McGill University (2013 - present) *
+!* *
+!* -> created: 2011 *
+!* -> latest changes: 2019/07/29 *
+!* *
+!* :: License :: *
+!* This program is free software: you can redistribute it and/or modify it under the *
+!* terms of the GNU General Public License as published by the Free Software *
+!* Foundation, either version 3 of the License, or (at your option) any later *
+!* version. *
+!* The AIOMFAC model code 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 General Public License for more *
+!* details. *
+!* You should have received a copy of the GNU General Public License along with this *
+!* program. If not, see . *
+!* *
+!****************************************************************************************
+SUBROUTINE ReadInputFile(filepath, filename, filepathout, ninpmax, maxpoints, unito, verbose, ncp, npoints, &
+ & warningind, errorind, filevalid, cpnameinp, cpsubg, T_K, composition, xinputtype)
+
+USE ModSystemProp, ONLY : topsubno
+
+IMPLICIT NONE
+
+!interface:
+CHARACTER(LEN=3000),INTENT(INOUT) :: filepath, filename, filepathout
+INTEGER(4),INTENT(IN) :: ninpmax, maxpoints
+INTEGER(4),INTENT(INOUT) :: unito
+LOGICAL(4),INTENT(IN) :: verbose
+INTEGER(4),INTENT(OUT) :: ncp, npoints
+INTEGER(4),INTENT(INOUT) :: warningind, errorind
+LOGICAL(4),INTENT(OUT) :: filevalid
+CHARACTER(LEN=200),DIMENSION(:),INTENT(OUT) :: cpnameinp !list of assigned component names (from input file)
+INTEGER(4),DIMENSION(:,:),INTENT(OUT) :: cpsubg !list of input component subgroups and corresponding subgroup quantities
+REAL(8),DIMENSION(:),INTENT(OUT) :: T_K !temperature of data points in Kelvin
+REAL(8),DIMENSION(:,:),INTENT(OUT) :: composition !array of mixture composition points for which calculations should be run
+LOGICAL(4),INTENT(OUT) :: xinputtype
+!--
+!local variables:
+CHARACTER(LEN=:),ALLOCATABLE :: cn !this assumes a maximum four-digit component number in the system (max. 9999); to be adjusted otherwise.
+CHARACTER(LEN=4) :: dashes, equalsigns, pluses
+CHARACTER(LEN=20) :: dummy, cnformat
+CHARACTER(LEN=50) :: txtcheck, inpfolder
+CHARACTER(LEN=3000) :: errlogfile, fname
+CHARACTER(LEN=20),DIMENSION(ninpmax) :: txtarray
+INTEGER(4) :: cpno, i, inpfilesize, istat, k, kinpf, qty, subg, unitx
+LOGICAL(4) :: fileopened, fileexists
+!...................................................................................
+
+!initialize variables
+ncp = 0
+cpsubg = 0
+cpnameinp = "none"
+composition = 0.0D0
+T_K = 298.15D0 !default temperature in [K]
+dashes = "----"
+equalsigns = "===="
+pluses = "++++"
+fileexists = .false.
+filevalid = .false.
+
+!==== READ INPUT data ===========================================================
+
+!extract filepath from the input filepath (which could include a path to the directory):
+k = LEN_TRIM(filepath)
+IF (k < 1) THEN !no file was provided
+ WRITE(*,*) ""
+ WRITE(*,*) "ERROR from AIOMFAC-web: no input file was provided!"
+ WRITE(*,*) "An input file path needs to be provided via command line argument 2"
+ WRITE(*,*) ""
+ READ(*,*)
+ STOP
+ENDIF
+
+inpfolder = "Inputfiles" !"Inputfiles" in AIOMFAC-web
+kinpf = LEN_TRIM(inpfolder)
+i = INDEX(filepath, TRIM(inpfolder)//"/input")
+IF (i > 0 .AND. i < k) THEN !found a direcory path in UNIX file path style
+ filename = filepath(i+kinpf+1:k)
+ filepath = filepath(1:i+kinpf)
+ELSE
+ i = INDEX(filepath, TRIM(inpfolder)//"\input")
+ IF (i > 0 .AND. i < k) THEN !found a direcory path in WINDOWS file path style
+ filename = filepath(i+kinpf+1:k)
+ filepath = filepath(1:i+kinpf)
+ ELSE
+ filename = filepath(i+kinpf+1:k)
+ filepath = ""
+ ENDIF
+ENDIF
+!-----
+!check / create for associated output directory "Outputfiles":
+i = LEN_TRIM(filepath)
+filepathout = TRIM(filepath(1:i-(kinpf+1)))//"Outputfiles/"
+
+!use the name of the input file to create a corresponding output file name:
+i = INDEX(filename, ".txt") !returns starting position of string input within string filename
+!create an error-logfile associated with the input file name:
+errlogfile = "Errorlog_"//filename(i-4:)
+fname = TRIM(filepathout)//TRIM(errlogfile)
+OPEN (NEWUNIT = unito, FILE = fname, STATUS ='UNKNOWN') !unito is the error / logfile unit
+!-----
+!check if file exists and read its content if true:
+fname = TRIM(filepath)//TRIM(filename)
+INQUIRE(FILE = fname, EXIST = fileexists, SIZE = inpfilesize) !inpfilesize is the file size in [bytes]
+!Delete very large files that can only mean uploaded spam content and not actual input:
+IF (fileexists) THEN
+ IF (FLOAT(inpfilesize) > 50*(ninpmax +ninpmax*maxpoints) ) THEN !cannot be a valid input file
+ fname = TRIM(filepath)//TRIM(filename)
+ OPEN (NEWUNIT = unitx, FILE = fname, STATUS='OLD')
+ CLOSE(unitx, STATUS='DELETE') !close and delete the file
+ fileexists = .false.
+ ENDIF
+ENDIF
+
+!attempt to read a supposedly valid input file:
+IF (fileexists) THEN
+ IF (verbose) THEN
+ WRITE(unito,*) ""
+ WRITE(unito,*) "MESSAGE from AIOMFAC: input file exists."
+ ENDIF
+ fname = TRIM(filepath)//TRIM(filename)
+ OPEN (NEWUNIT = unitx, FILE = fname, IOSTAT=istat, ACTION='READ', STATUS='OLD')
+ IF (istat /= 0) THEN ! an error occurred
+ WRITE(unito,*) ""
+ WRITE(unito,*) "MESSAGE from AIOMFAC: an error occurred while trying to open the file! IOSTAT: ", istat
+ ENDIF
+ !validate file as a correct input text-file (no spam):
+ READ(unitx,*) dummy, dummy, dummy, txtcheck
+ IF (txtcheck(1:11) == "AIOMFAC-web") THEN
+ filevalid = .true.
+ IF (verbose) THEN
+ WRITE(unito,*) ""
+ WRITE(unito,*) "MESSAGE from AIOMFAC: input file has passed the first line text validation and will be read."
+ ENDIF
+ !read input data from file:
+ BACKSPACE unitx !jump back to beginning of record (to the beginning of the line)
+ READ(unitx,*) dummy, dummy, dummy, dummy, dummy !read first line
+ READ(unitx,*) !read empty line 2
+ READ(unitx,*) dummy, dummy !read line 3
+ READ(unitx,*) dummy !read line 4
+ ELSE !file not valid. It will be closed and deleted below
+ filevalid = .false.
+ IF (verbose) THEN
+ WRITE(unito,*) ""
+ WRITE(unito,*) "MESSAGE from AIOMFAC: input file does not pass first line text validation and will be deleted."
+ ENDIF
+ ENDIF
+ !loop over mixture components with variable numbers of subgroups to read (using inner loop):
+ IF (filevalid) THEN
+ IF (verbose) THEN
+ WRITE(unito,*) ""
+ WRITE(unito,*) "MESSAGE from AIOMFAC: reading component data from input file."
+ ENDIF
+ k = MAX(2, CEILING(LOG10(REAL(ninpmax))) ) !determine order of magnitude digits
+ WRITE(dummy,'(I0)') k
+ cnformat = "(I"//TRIM(dummy)//"."//TRIM(dummy)//")" !variable format specifier
+ ALLOCATE(CHARACTER(LEN=k) :: cn)
+ DO ncp = 1,ninpmax
+ IF (verbose .AND. istat /= 0) THEN
+ WRITE(unito,*) ""
+ WRITE(unito,*) "MESSAGE from AIOMFAC: file end found in input file at ncp = ", ncp
+ ENDIF
+ READ(unitx,*,IOSTAT=istat) txtcheck !read only first argument on this line for subsequent check
+ IF (txtcheck(1:4) == pluses .OR. istat /= 0) THEN !"++++" indicates end of this components definition part
+ EXIT !exit ncp DO-loop
+ ELSE !in this case, argument 3 of txtcheck is the component no.:
+ BACKSPACE unitx !jump back to beginning of record (to the beginning of the line)
+ READ(unitx,*) dummy, dummy, cpno !read line with component's number
+ ENDIF
+ READ(unitx,*) dummy, dummy, cpnameinp(cpno) !read line with component's name
+ IF (LEN_TRIM(cpnameinp(cpno)) < 1) THEN !no component name was assigned, generate a default name
+ WRITE(cn,cnformat) cpno
+ cpnameinp(cpno) = "cp_"//cn
+ ENDIF
+ IF (verbose) THEN
+ WRITE(unito,*) ""
+ WRITE(unito,*) "MESSAGE from AIOMFAC: found a component cpno =", cpno
+ ENDIF
+ DO !until exit
+ READ(unitx,*,IOSTAT=istat) txtcheck !read only first argument on this line
+ !check whether another subgroup is present or not
+ IF (txtcheck(1:4) == dashes .OR. istat /= 0) THEN !"----" indicates no more subgroups of this component
+ EXIT !leave the inner DO-loop
+ ELSE !else continue reading the next subgroup
+ BACKSPACE unitx
+ READ(unitx,*) dummy, dummy, dummy, subg, qty !continue reading line with subgroup no. and corresp. quantity
+ cpsubg(cpno,subg) = cpsubg(cpno,subg)+qty
+ ENDIF
+ ENDDO
+ ENDDO
+ ncp = ncp-1 !ncp-1 is the number of different components in this mixture
+ IF (verbose) THEN
+ WRITE(unito,*) ""
+ WRITE(unito,*) "MESSAGE from AIOMFAC: component data read."
+ WRITE(unito,*) "MESSAGE from AIOMFAC: number of components: ", ncp
+ ENDIF
+ IF (ncp == ninpmax) THEN
+ READ(unitx,*,IOSTAT=istat) txtcheck !read only first argument on this line for subsequent check
+ IF (txtcheck(1:4) /= pluses) THEN
+ errorind = 34
+ filevalid = .false.
+ WRITE(*,*) "AIOMFAC ERROR 34: maximum number of input components reached while reading input file."
+ WRITE(*,*) "AIOMFAC ERROR 34: check whether ninpmax value is too small."
+ WRITE(unito,*) ""
+ WRITE(unito,*) "AIOMFAC ERROR 34: maximum number of input components reached while reading input file."
+ ENDIF
+ ENDIF
+ IF (filevalid) THEN
+ READ(unitx,*) dummy, dummy !read mixture composition line
+ READ(unitx,*) dummy, dummy, i !read mass fraction? line
+ IF (i == 1) THEN !composition in mass fractions
+ xinputtype = .false.
+ ELSE !composition in mole fractions
+ xinputtype = .true.
+ ENDIF
+ READ(unitx,*) dummy, dummy, i !read mole fraction? line
+ READ(unitx,*) dummy !read dashes line
+ !now read the lines with the mixture composition points:
+ READ(unitx,*) txtarray(1:ncp) !read header line of composition table
+ DO npoints = 1,maxpoints !or until exit
+ READ(unitx,*,IOSTAT=istat) txtcheck !read only the first argument on this line
+ IF (txtcheck(1:4) == equalsigns .OR. istat /= 0) THEN !"====" indicates no more composition points (and last line of input file)
+ EXIT !leave the DO-loop (normal exit point)
+ ELSE IF (IACHAR(txtcheck(1:1)) > 47 .AND. IACHAR(txtcheck(1:1)) < 58) THEN !validate that the data is actual intended input and not some sort of text field spam).
+ BACKSPACE unitx
+ READ(unitx,*) txtcheck, dummy !read only the first two arguments on this line
+ IF (IACHAR(dummy(1:1)) > 47 .AND. IACHAR(dummy(1:1)) < 58) THEN !it is a number
+ BACKSPACE unitx
+ READ(unitx,*) i, T_K(npoints), composition(npoints,2:ncp) !read the temperature in [K] and composition values of the components [2:ncp] into the array
+ ELSE
+ IF (npoints == 1) THEN
+ filevalid = .false. !file is not valid due to incorrect input in text field
+ EXIT
+ ELSE
+ warningind = 31
+ EXIT !abnormal exit point of loop due to non-number entries at a certain line in the text field
+ ENDIF
+ ENDIF
+ ELSE
+ IF (npoints == 1) THEN
+ filevalid = .false. !file is not valid due to incorrect input in text field
+ EXIT
+ ELSE
+ warningind = 31
+ EXIT !abnormal exit point of loop due to non-number entries at a certain line in the text field
+ ENDIF
+ ENDIF
+ ENDDO
+ npoints = npoints-1 !the number of composition points
+ IF (verbose) THEN
+ WRITE(unito,*) ""
+ WRITE(unito,*) "MESSAGE from AIOMFAC: composition points read."
+ WRITE(unito,*) "MESSAGE from AIOMFAC: number of points: ", npoints
+ ENDIF
+ IF (.NOT. filevalid) THEN
+ !close and delete file from server:
+ CLOSE(unitx, STATUS = 'DELETE')
+ ELSE
+ CLOSE(unitx)
+ ENDIF
+ !assign component 01 its composition values based on the rest of the component's data:
+ DO i = 1,npoints
+ composition(i,1) = 1.0D0-SUM(composition(i,2:ncp))
+ ENDDO
+ ENDIF !filevalid2
+ ENDIF !filevalid1
+ELSE
+ WRITE(unito,*) ""
+ WRITE(unito,*) "ERROR in AIOMFAC: Input file does not exist at expected location: ", TRIM(filename)
+ WRITE(unito,*) ""
+ENDIF !fileexists
+
+IF (fileexists) THEN
+ IF (.NOT. filevalid) THEN !input file contains errors or is completely invalid (submitted spam etc.)
+ IF (verbose) THEN
+ WRITE(unito,*) ""
+ WRITE(unito,*) "MESSAGE from AIOMFAC: input file did not pass full validation and may be a spam file."
+ WRITE(unito,*) "MESSAGE from AIOMFAC: the input file will be deleted to prevent spam files and malicious code from occupying the server."
+ WRITE(unito,*) ""
+ ENDIF
+ INQUIRE(FILE = fname, OPENED = fileopened)
+ IF (errorind == 0) THEN
+ errorind = 32
+ !close and delete file from server:
+ IF (fileopened) THEN
+ CLOSE(unitx, STATUS = 'DELETE')
+ ENDIF
+ ELSE
+ IF (fileopened) THEN
+ CLOSE(unitx)
+ ENDIF
+ ENDIF
+ ELSE
+ !check for any valid points for calculations before initializing AIOMFAC:
+ IF (npoints < 1) THEN !no points input for calculations
+ errorind = 33 !no input in text field..
+ IF (verbose) THEN
+ WRITE(unito,*) ""
+ WRITE(unito,*) "MESSAGE from AIOMFAC: no composition points have been entered in the text field. There is nothing to calculate."
+ WRITE(unito,*) ""
+ ENDIF
+ filevalid = .false.
+ ENDIF
+ ENDIF
+ENDIF
+
+END SUBROUTINE ReadInputFile
\ No newline at end of file
diff --git a/FortranCode/RepErrorWarning.f90 b/FortranCode/RepErrorWarning.f90
index 6016cdc..94b0842 100644
--- a/FortranCode/RepErrorWarning.f90
+++ b/FortranCode/RepErrorWarning.f90
@@ -9,8 +9,8 @@
!* *
!* -> created: 2011 *
!* -> latest changes: 2018/08/08 *
-!* *
-!* :: License :: *
+!* *
+!* :: License :: *
!* This program is free software: you can redistribute it and/or modify it under the *
!* terms of the GNU General Public License as published by the Free Software *
!* Foundation, either version 3 of the License, or (at your option) any later *
@@ -127,12 +127,23 @@ SUBROUTINE RepErrorWarning(unito, errorflagmix, warningflag, errorflagcalc, poin
WRITE(unito,*) "AIOMFAC WARNING 11: Temperature range related."
WRITE(unito,*) "At least one data point has a set temperature outside of"
WRITE(unito,*) "the recommended range for model calculations of "
- WRITE(unito,*) "electrolyte-free organic mixtures. This may be intended, "
+ WRITE(unito,*) "electrolyte-free organic mixtures. This may be intended,"
WRITE(unito,*) "but caution is advised as AIOMFAC is not designed to "
WRITE(unito,*) "perform well at this temperature."
WRITE(unito,*) "Data point no.: ", pointi
WRITE(unito,*) "======================================================="
WRITE(unito,*) ""
+ CASE(16)
+ WRITE(unito,*) ""
+ WRITE(unito,*) "======================================================="
+ WRITE(unito,*) "AIOMFAC-VISC WARNING 16: Mixture viscosity issue. "
+ WRITE(unito,*) "Note that mixture viscosity is currently not computed "
+ WRITE(unito,*) "for electrolyte-containing mixtures. Therefore, an "
+ WRITE(unito,*) "unrealistic mixture viscosity of "
+ WRITE(unito,*) "log_10(eta/[Pa.s]) = -999.999 is output. "
+ WRITE(unito,*) "Data point no.: ", pointi
+ WRITE(unito,*) "======================================================="
+ WRITE(unito,*) ""
CASE DEFAULT
WRITE(unito,*) ""
WRITE(unito,*) "======================================================="
diff --git a/FortranCode/SubModDefSystem.f90 b/FortranCode/SubModDefSystem.f90
index 4e57336..7488df3 100644
--- a/FortranCode/SubModDefSystem.f90
+++ b/FortranCode/SubModDefSystem.f90
@@ -12,8 +12,8 @@
!* *
!* -> created: 2005 *
!* -> latest changes: 2018/05/31 *
-!* *
-!* :: License :: *
+!* *
+!* :: License :: *
!* This program is free software: you can redistribute it and/or modify it under the *
!* terms of the GNU General Public License as published by the Free Software *
!* Foundation, either version 3 of the License, or (at your option) any later *
@@ -66,7 +66,7 @@ MODULE SUBROUTINE SetSystem(ndi, datafromfile, ninp, cpnameinp, cpsubginp)
LOGICAL(4),INTENT(IN) :: datafromfile !set .true. if the system components are provided from an input file
INTEGER(4),INTENT(IN) :: ninp !value known at input only in the case of input from a file
!optional input arguments:
- CHARACTER(LEN=60),DIMENSION(ninp),INTENT(IN), OPTIONAL :: cpnameinp
+ CHARACTER(LEN=200),DIMENSION(ninp),INTENT(IN), OPTIONAL :: cpnameinp
INTEGER(4),DIMENSION(ninp,topsubno),INTENT(IN), OPTIONAL :: cpsubginp
!...
!local variables
@@ -233,7 +233,6 @@ MODULE SUBROUTINE definemixtures(ndi, ninputcomp, compID, cpsubg)
INTEGER(4),DIMENSION(200) :: SolvSubs2 !temporary array for solvent subgroups (we allow a maximum of 200)
INTEGER(4),DIMENSION(ninputcomp*2) :: ElectSubs2
INTEGER(4),DIMENSION(201:topsubno) :: ElectPos
- CHARACTER(LEN=3) :: cn
LOGICAL(4) :: already1, already2
!--------------------------------------