Greetings everyone

I am working on a modeling project that involves some complex computations. I transitioned a model from R to Fortran due to performance concerns that hindered the efficiency of the original model being to slow in R. Within this Fortran environment, I am in need of robust root-finding algorithms. While exploring, I came across the use of root-fortran library for finding root in modern Fortran. However on translating the r code to fortran and solving the code using similar tolerance as used in R code I have not been getting the same values as in the R code. The R function is a bit complex and translated to fortran as below

function Conc_calc(zHV) result(result)

real(Kind=8), intent(in) :: zHV

real(Kind=8) :: A, B, Kw, K1, K2, K3, H1, H2, zH

real(Kind=8) :: zH0, pH0, pHI, buff_L

real(kind=8), dimension(12) :: vec1, Vec2, vec3, vec4, vec5, vec6, vec7

real(kind=8), dimension(12) :: vec8, vec9, vec10, vec11, vec12

real(kind=8) :: Amatrix(12,12), Bvector(12)

real(kind=8) :: roots(12)

real(kind=8) :: result`integer :: n=12 ! specify the size of the matrix A = x10 + x20 + x30 + InN - OutN B = y10 + y20 + y30 + y40 + InC Kw = 10.0**(-14.0) ! Kw = 10.0**(-14.0) at 25 C !Print *, 'zHV', Tf, ti K1 = 5.67 * (10.0**(-10.0)) * exp(-6286.0 * ((1.0 / Tf) - (1.0 / 298.15))) !defines the first H2CO3 dissociation into HCO3- K2 = (10.0**(-1.0 * (((2902.39 / Tf) + (0.02379 * Tf) - 6.4980)))) !defines the second HCO3- dissociation into CO3-- K3 = (10.0**(-1.0 * (((3404.71 / Tf) + (0.032786 * Tf) - 14.8435)))) !defines the third H2CO3-- dissociation into H2O and co2 H1 = (56.0 * exp(4092.0 * ((1.0 / Tf) - (1.0 / 298.15)))) / ((1.0 / 8.314) * 0.001 * (1.0 / Tf) * & (1.013 * 10.0 ** 5)) H2 = ((1.0 / 29.41) * exp(2400.0 * ((1.0 / Tf) - (1.0 / 298.15)))) / ((1.0 / 8.314) * 0.001 * & (1.0 / Tf) * (1.013 * 10.0 ** 5)) zH = zHV * V if(ti == 2 .and. ti3 == 3) then write(*,*) 'zH is ',zH, ti, ti3 end if vec1 = (/0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 1.0_wp, 0.0_wp, -zH, 0.0_wp, & 0.0_wp, (1.0_wp + H1 * (V / Vg)), 0.0_wp/) ! x1 vec2 = (/0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 1.0_wp, 0.0_wp, 0.0_wp, K1* V, 0.0_wp, 0.0_wp, H1 * (V / Vg), 0.0_wp/) ! x2 vec6 = (/0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 1.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp/) vec3 = (/0.0_wp, 0.0_wp, 1.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, -zH, 0.0_wp, 0.0_wp, H2* (V / Vg)/) ! y1 vec4 = (/0.0_wp, 1.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, K2*V, -zH, 0.0_wp, H2 * (V / Vg)/) ! y2 vec5 = (/1.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, & K3* V, 0.0_wp, (1.0_wp + H2 * (V / Vg))/) ! y3 vec7 = (/0.0_wp, 0.0_wp, 0.0_wp, 1.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp/) ! y4 vec8 = (/1.0_wp, 0.0_wp, 0.0_wp, -1.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp/) ! RH2 vec9 = (/-1.0_wp, 1.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp/) ! RK3 vec10 = (/0.0_wp, 1.0_wp, -1.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp/) ! RK2 vec11 = (/0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 1.0_wp, -1.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp/) ! RK1 vec12 = (/0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 1.0_wp, -1.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp/) ! RH1 Amatrix = reshape([vec1, vec2, vec6, vec3, vec4, vec5, vec7, vec8, vec9, vec10, & vec11, vec12], shape(Amatrix)) ! Set the column names !character(2) :: colnames(12) = ["x1", "x2", "x3", "y1", "y2", "y3", "y4", "RH2", "RK3", & ! "RK2", "RK1", "RH1"] Bvector = (/y30, y20 + InC, y10, y40 - OutC, x20 + InN, x10, x30 - OutN, 0.0_wp, 0.0_wp, 0.0_wp, & H1 * A * (V / Vg), H2 * B * (V / Vg)/) ! const call solve_system(n, Amatrix, Bvector, roots) !I called lapack to solve the system zH0 = zHV0 * V pH0 = -log10(zH0 / V) pHI = -log10(zH / V) !I renamed the line pH to PHI result =(zH0 - InC + (-roots(9) + roots(10) + roots(11)) + (buff_L * (pHI - pH0))) - ((zH)) end function Conc_calc`

I used the root module of root fortran and I called it as below

call root_scalar(âbrentâ,Conc_calc, 0.0_wp, 1.0_wp, zHV, result, iflag, atol=1.0_wp,Ftol = 0.10_wp**(0.25_wp))

However for the last 3 days I have been stuck as I havent been getting the same result as the R code. I was wondering if I got a step wrong or not so I started to print every single variable. I noticed that when I print zH which is within the function I am finding the root of in both fortran and R I get the value to be in fortran as below

.\my_program.exe zH is0.0000000000000000

zH is 57.641325192598238

zH is 28.820662596299119

zH is 14.410331298149559

zH is 7.2051656490747797

zH is 3.6025828245373899

zH is 1.8012914122686949

zH is 0.90064570613434747

zH is 0.45032285306717373

zH is 0.22516142653358687

zH is 0.11258071326679343

zH is 5.6290356633396717E-002

zH is 2.8145178316698358E-002

zH is 1.4072589158349179E-002

zH is 7.0362945791745896E-003

zH is 3.5181472895872948E-003

zH is 1.7590736447936474E-003the last line in fortran

where as in the R code I got values lower than this

1] zH is 0"

[1] " zH is 57.6413255382803"

[1] " zH is 28.8206627691402"

[1] " zH is 14.4103313845701"

[1] " zH is 7.20516569228504"

[1] " zH is 3.60258284614252"

[1] " zH is 1.80129142307126"

[1] " zH is 0.90064571153563"

[1] " zH is 0.450322855767815"

[1] " zH is 0.225161427883908"

[1] " zH is 0.112580713941954"

[1] " zH is 0.0562903569709769"

[1] " zH is 0.0281451784854884"

[1] " zH is 0.0140725892427442"

[1] " zH is 0.00703629462137211"

[1] " zH is 0.00351814731068605"

[1]" zH is 0.00175907365534303"where the fortran ended

[1] " zH is 0.000879536827671514"

[1] " zH is 0.000489560168170512"

[1] " zH is 0.000674886029990813"

[1] " zH is 0.000650298707008735"

[1] " zH is 0.00064673969984793"

[1] " zH is 0.000646746327972608"

[1] " zH is 0.000646746309843868"

[1] " zH is 0.000646746309843775"

[1] " zH is 0.000646746309843776"

[1] " zH is 0.000646746309843776"

After using the one specified in the R code since in R, the default setting is. tol = .Machine$double.eps^0.25.

I have really been tweaking the ftol value which has been changing the values but stopped at the value above in the call rootscalar**(used 1.0_wp**(0.25_wp)** )in fortran. I am working with ecological data and this step is an important one to finish up with the translation.

I would be very grateful if I get suggestions on the next step to take to figure out this as I have been stuck on this step for almost a week. Many thanks in advance