1 | !> \file mocsy_phsolvers.f90 |
---|
2 | !! \BRIEF |
---|
3 | !> Module with routines needed to solve pH-total alkalinity equation (Munhoven, 2013, GMD) |
---|
4 | MODULE mocsy_phsolvers |
---|
5 | ! Module of fastest solvers from Munhoven (2013, Geosci. Model Dev., 6, 1367-1388) |
---|
6 | ! ! Taken from SolveSAPHE (mod_phsolvers.F90) & adapted very slightly for use with mocsy |
---|
7 | ! ! SolveSaphe is distributed under the GNU Lesser General Public License |
---|
8 | ! ! mocsy is distributed under the MIT License |
---|
9 | ! |
---|
10 | ! Modifications J. C. Orr, LSCE/IPSL, CEA-CNRS-UVSQ, France, 11 Sep 2014: |
---|
11 | ! 1) kept only the 3 fastest solvers (atgen, atsec, atfast) and routines which they call |
---|
12 | ! 2) reduced vertical white space: deleted many blank lines & comment lines that served as divisions |
---|
13 | ! 3) converted name from .F90 to .f90, deleting a few optional preprocesse if statements |
---|
14 | ! 4) read in mocsy computed equilibrium constants (as arguments) instead of USE MOD_CHEMCONST |
---|
15 | ! 5) converted routine names from upper case to lower case |
---|
16 | ! 6) commented out arguments and equations for NH4 and H2S acid systems |
---|
17 | |
---|
18 | USE mocsy_singledouble |
---|
19 | IMPLICIT NONE |
---|
20 | |
---|
21 | ! General parameters |
---|
22 | REAL(KIND=wp), PARAMETER :: pp_rdel_ah_target = 1.E-8_wp |
---|
23 | REAL(KIND=wp), PARAMETER :: pp_ln10 = 2.302585092994045684018_wp |
---|
24 | |
---|
25 | ! Maximum number of iterations for each method |
---|
26 | INTEGER, PARAMETER :: jp_maxniter_atgen = 50 |
---|
27 | INTEGER, PARAMETER :: jp_maxniter_atsec = 50 |
---|
28 | INTEGER, PARAMETER :: jp_maxniter_atfast = 50 |
---|
29 | |
---|
30 | ! Bookkeeping variables for each method |
---|
31 | ! - SOLVE_AT_GENERAL |
---|
32 | INTEGER :: niter_atgen = jp_maxniter_atgen |
---|
33 | |
---|
34 | ! - SOLVE_AT_GENERAL_SEC |
---|
35 | INTEGER :: niter_atsec = jp_maxniter_atsec |
---|
36 | |
---|
37 | ! - SOLVE_AT_FAST (variant of SOLVE_AT_GENERAL w/o bracketing |
---|
38 | INTEGER :: niter_atfast = jp_maxniter_atfast |
---|
39 | |
---|
40 | ! Keep the following functions private to avoid conflicts with |
---|
41 | ! other modules that provide similar ones. |
---|
42 | !PRIVATE AHINI_FOR_AT |
---|
43 | |
---|
44 | CONTAINS |
---|
45 | !=============================================================================== |
---|
46 | SUBROUTINE anw_infsup(p_dictot, p_bortot, & |
---|
47 | p_po4tot, p_siltot, & |
---|
48 | p_so4tot, p_flutot, & |
---|
49 | p_alknw_inf, p_alknw_sup) |
---|
50 | |
---|
51 | ! Subroutine returns the lower and upper bounds of "non-water-selfionization" |
---|
52 | ! contributions to total alkalinity (the infimum and the supremum), i.e |
---|
53 | ! inf(TA - [OH-] + [H+]) and sup(TA - [OH-] + [H+]) |
---|
54 | |
---|
55 | USE mocsy_singledouble |
---|
56 | IMPLICIT NONE |
---|
57 | |
---|
58 | ! Argument variables |
---|
59 | REAL(KIND=wp), INTENT(IN) :: p_dictot |
---|
60 | REAL(KIND=wp), INTENT(IN) :: p_bortot |
---|
61 | REAL(KIND=wp), INTENT(IN) :: p_po4tot |
---|
62 | REAL(KIND=wp), INTENT(IN) :: p_siltot |
---|
63 | !REAL(KIND=wp), INTENT(IN) :: p_nh4tot |
---|
64 | !REAL(KIND=wp), INTENT(IN) :: p_h2stot |
---|
65 | REAL(KIND=wp), INTENT(IN) :: p_so4tot |
---|
66 | REAL(KIND=wp), INTENT(IN) :: p_flutot |
---|
67 | REAL(KIND=wp), INTENT(OUT) :: p_alknw_inf |
---|
68 | REAL(KIND=wp), INTENT(OUT) :: p_alknw_sup |
---|
69 | |
---|
70 | p_alknw_inf = -p_po4tot - p_so4tot - p_flutot |
---|
71 | p_alknw_sup = p_dictot + p_dictot + p_bortot & |
---|
72 | + p_po4tot + p_po4tot + p_siltot !& |
---|
73 | ! + p_nh4tot + p_h2stot |
---|
74 | |
---|
75 | RETURN |
---|
76 | END SUBROUTINE anw_infsup |
---|
77 | |
---|
78 | !=============================================================================== |
---|
79 | |
---|
80 | FUNCTION equation_at(p_alktot, p_h, p_dictot, p_bortot, & |
---|
81 | p_po4tot, p_siltot, & |
---|
82 | p_so4tot, p_flutot, & |
---|
83 | K0, K1, K2, Kb, Kw, Ks, Kf, K1p, K2p, K3p, Ksi, & |
---|
84 | p_deriveqn) |
---|
85 | |
---|
86 | USE mocsy_singledouble |
---|
87 | IMPLICIT NONE |
---|
88 | REAL(KIND=wp) :: equation_at |
---|
89 | |
---|
90 | ! Argument variables |
---|
91 | REAL(KIND=wp), INTENT(IN) :: p_alktot |
---|
92 | REAL(KIND=wp), INTENT(IN) :: p_h |
---|
93 | REAL(KIND=wp), INTENT(IN) :: p_dictot |
---|
94 | REAL(KIND=wp), INTENT(IN) :: p_bortot |
---|
95 | REAL(KIND=wp), INTENT(IN) :: p_po4tot |
---|
96 | REAL(KIND=wp), INTENT(IN) :: p_siltot |
---|
97 | !REAL(KIND=wp), INTENT(IN) :: p_nh4tot |
---|
98 | !REAL(KIND=wp), INTENT(IN) :: p_h2stot |
---|
99 | REAL(KIND=wp), INTENT(IN) :: p_so4tot |
---|
100 | REAL(KIND=wp), INTENT(IN) :: p_flutot |
---|
101 | REAL(KIND=wp), INTENT(IN) :: K0, K1, K2, Kb, Kw, Ks, Kf |
---|
102 | REAL(KIND=wp), INTENT(IN) :: K1p, K2p, K3p, Ksi |
---|
103 | REAL(KIND=wp), INTENT(OUT), OPTIONAL :: p_deriveqn |
---|
104 | |
---|
105 | ! Local variables |
---|
106 | !----------------- |
---|
107 | REAL(KIND=wp) :: znumer_dic, zdnumer_dic, zdenom_dic, zalk_dic, zdalk_dic |
---|
108 | REAL(KIND=wp) :: znumer_bor, zdnumer_bor, zdenom_bor, zalk_bor, zdalk_bor |
---|
109 | REAL(KIND=wp) :: znumer_po4, zdnumer_po4, zdenom_po4, zalk_po4, zdalk_po4 |
---|
110 | REAL(KIND=wp) :: znumer_sil, zdnumer_sil, zdenom_sil, zalk_sil, zdalk_sil |
---|
111 | REAL(KIND=wp) :: znumer_nh4, zdnumer_nh4, zdenom_nh4, zalk_nh4, zdalk_nh4 |
---|
112 | REAL(KIND=wp) :: znumer_h2s, zdnumer_h2s, zdenom_h2s, zalk_h2s, zdalk_h2s |
---|
113 | REAL(KIND=wp) :: znumer_so4, zdnumer_so4, zdenom_so4, zalk_so4, zdalk_so4 |
---|
114 | REAL(KIND=wp) :: znumer_flu, zdnumer_flu, zdenom_flu, zalk_flu, zdalk_flu |
---|
115 | REAL(KIND=wp) :: zalk_wat, zdalk_wat |
---|
116 | REAL(KIND=wp) :: aphscale |
---|
117 | |
---|
118 | ! TOTAL H+ scale: conversion factor for Htot = aphscale * Hfree |
---|
119 | aphscale = 1._wp + p_so4tot/Ks |
---|
120 | |
---|
121 | ! H2CO3 - HCO3 - CO3 : n=2, m=0 |
---|
122 | znumer_dic = 2._wp*K1*K2 + p_h* K1 |
---|
123 | zdenom_dic = K1*K2 + p_h*( K1 + p_h) |
---|
124 | zalk_dic = p_dictot * (znumer_dic/zdenom_dic) |
---|
125 | |
---|
126 | ! B(OH)3 - B(OH)4 : n=1, m=0 |
---|
127 | znumer_bor = Kb |
---|
128 | zdenom_bor = Kb + p_h |
---|
129 | zalk_bor = p_bortot * (znumer_bor/zdenom_bor) |
---|
130 | |
---|
131 | ! H3PO4 - H2PO4 - HPO4 - PO4 : n=3, m=1 |
---|
132 | znumer_po4 = 3._wp*K1p*K2p*K3p + p_h*(2._wp*K1p*K2p + p_h* K1p) |
---|
133 | zdenom_po4 = K1p*K2p*K3p + p_h*( K1p*K2p + p_h*(K1p + p_h)) |
---|
134 | zalk_po4 = p_po4tot * (znumer_po4/zdenom_po4 - 1._wp) ! Zero level of H3PO4 = 1 |
---|
135 | |
---|
136 | ! H4SiO4 - H3SiO4 : n=1, m=0 |
---|
137 | znumer_sil = Ksi |
---|
138 | zdenom_sil = Ksi + p_h |
---|
139 | zalk_sil = p_siltot * (znumer_sil/zdenom_sil) |
---|
140 | |
---|
141 | ! NH4 - NH3 : n=1, m=0 |
---|
142 | !znumer_nh4 = api1_nh4 |
---|
143 | !zdenom_nh4 = api1_nh4 + p_h |
---|
144 | !zalk_nh4 = p_nh4tot * (znumer_nh4/zdenom_nh4) |
---|
145 | ! Note: api1_nh4 = Knh4 |
---|
146 | |
---|
147 | ! H2S - HS : n=1, m=0 |
---|
148 | !znumer_h2s = api1_h2s |
---|
149 | !zdenom_h2s = api1_h2s + p_h |
---|
150 | !zalk_h2s = p_h2stot * (znumer_h2s/zdenom_h2s) |
---|
151 | ! Note: api1_h2s = Kh2s |
---|
152 | |
---|
153 | ! HSO4 - SO4 : n=1, m=1 |
---|
154 | znumer_so4 = Ks |
---|
155 | zdenom_so4 = Ks + p_h |
---|
156 | zalk_so4 = p_so4tot * (znumer_so4/zdenom_so4 - 1._wp) |
---|
157 | |
---|
158 | ! HF - F : n=1, m=1 |
---|
159 | znumer_flu = Kf |
---|
160 | zdenom_flu = Kf + p_h |
---|
161 | zalk_flu = p_flutot * (znumer_flu/zdenom_flu - 1._wp) |
---|
162 | |
---|
163 | ! H2O - OH |
---|
164 | zalk_wat = Kw/p_h - p_h/aphscale |
---|
165 | |
---|
166 | equation_at = zalk_dic + zalk_bor + zalk_po4 + zalk_sil & |
---|
167 | + zalk_so4 + zalk_flu & |
---|
168 | + zalk_wat - p_alktot |
---|
169 | |
---|
170 | IF(PRESENT(p_deriveqn)) THEN |
---|
171 | ! H2CO3 - HCO3 - CO3 : n=2 |
---|
172 | zdnumer_dic = K1*K1*K2 + p_h*(4._wp*K1*K2 & |
---|
173 | + p_h* K1 ) |
---|
174 | zdalk_dic = -p_dictot*(zdnumer_dic/zdenom_dic**2) |
---|
175 | |
---|
176 | ! B(OH)3 - B(OH)4 : n=1 |
---|
177 | zdnumer_bor = Kb |
---|
178 | zdalk_bor = -p_bortot*(zdnumer_bor/zdenom_bor**2) |
---|
179 | |
---|
180 | ! H3PO4 - H2PO4 - HPO4 - PO4 : n=3 |
---|
181 | zdnumer_po4 = K1p*K2p*K1p*K2p*K3p + p_h*(4._wp*K1p*K1p*K2p*K3p & |
---|
182 | + p_h*(9._wp*K1p*K2p*K3p + K1p*K1p*K2p & |
---|
183 | + p_h*(4._wp*K1p*K2p & |
---|
184 | + p_h* K1p))) |
---|
185 | zdalk_po4 = -p_po4tot * (zdnumer_po4/zdenom_po4**2) |
---|
186 | |
---|
187 | ! H4SiO4 - H3SiO4 : n=1 |
---|
188 | zdnumer_sil = Ksi |
---|
189 | zdalk_sil = -p_siltot * (zdnumer_sil/zdenom_sil**2) |
---|
190 | |
---|
191 | ! ! NH4 - NH3 : n=1 |
---|
192 | ! zdnumer_nh4 = Knh4 |
---|
193 | ! zdalk_nh4 = -p_nh4tot * (zdnumer_nh4/zdenom_nh4**2) |
---|
194 | |
---|
195 | ! ! H2S - HS : n=1 |
---|
196 | ! zdnumer_h2s = api1_h2s |
---|
197 | ! zdalk_h2s = -p_h2stot * (zdnumer_h2s/zdenom_h2s**2) |
---|
198 | |
---|
199 | ! HSO4 - SO4 : n=1 |
---|
200 | zdnumer_so4 = Ks |
---|
201 | zdalk_so4 = -p_so4tot * (zdnumer_so4/zdenom_so4**2) |
---|
202 | |
---|
203 | ! HF - F : n=1 |
---|
204 | zdnumer_flu = Kf |
---|
205 | zdalk_flu = -p_flutot * (zdnumer_flu/zdenom_flu**2) |
---|
206 | |
---|
207 | ! p_deriveqn = zdalk_dic + zdalk_bor + zdalk_po4 + zdalk_sil & |
---|
208 | ! + zdalk_nh4 + zdalk_h2s + zdalk_so4 + zdalk_flu & |
---|
209 | ! - Kw/p_h**2 - 1._wp/aphscale |
---|
210 | p_deriveqn = zdalk_dic + zdalk_bor + zdalk_po4 + zdalk_sil & |
---|
211 | + zdalk_so4 + zdalk_flu & |
---|
212 | - Kw/p_h**2 - 1._wp/aphscale |
---|
213 | ENDIF |
---|
214 | RETURN |
---|
215 | END FUNCTION equation_at |
---|
216 | |
---|
217 | !=============================================================================== |
---|
218 | |
---|
219 | SUBROUTINE ahini_for_at(p_alkcb, p_dictot, p_bortot, K1, K2, Kb, p_hini) |
---|
220 | |
---|
221 | ! Subroutine returns the root for the 2nd order approximation of the |
---|
222 | ! DIC -- B_T -- A_CB equation for [H+] (reformulated as a cubic polynomial) |
---|
223 | ! around the local minimum, if it exists. |
---|
224 | |
---|
225 | ! Returns * 1E-03_wp if p_alkcb <= 0 |
---|
226 | ! * 1E-10_wp if p_alkcb >= 2*p_dictot + p_bortot |
---|
227 | ! * 1E-07_wp if 0 < p_alkcb < 2*p_dictot + p_bortot |
---|
228 | ! and the 2nd order approximation does not have a solution |
---|
229 | |
---|
230 | !USE MOD_CHEMCONST, ONLY : api1_dic, api2_dic, api1_bor |
---|
231 | |
---|
232 | USE mocsy_singledouble |
---|
233 | IMPLICIT NONE |
---|
234 | |
---|
235 | ! Argument variables |
---|
236 | !-------------------- |
---|
237 | REAL(KIND=wp), INTENT(IN) :: p_alkcb, p_dictot, p_bortot |
---|
238 | REAL(KIND=wp), INTENT(IN) :: K1, K2, Kb |
---|
239 | REAL(KIND=wp), INTENT(OUT) :: p_hini |
---|
240 | |
---|
241 | ! Local variables |
---|
242 | !----------------- |
---|
243 | REAL(KIND=wp) :: zca, zba |
---|
244 | REAL(KIND=wp) :: zd, zsqrtd, zhmin |
---|
245 | REAL(KIND=wp) :: za2, za1, za0 |
---|
246 | |
---|
247 | IF (p_alkcb <= 0._wp) THEN |
---|
248 | p_hini = 1.e-3_wp |
---|
249 | ELSEIF (p_alkcb >= (2._wp*p_dictot + p_bortot)) THEN |
---|
250 | p_hini = 1.e-10_wp |
---|
251 | ELSE |
---|
252 | zca = p_dictot/p_alkcb |
---|
253 | zba = p_bortot/p_alkcb |
---|
254 | |
---|
255 | ! Coefficients of the cubic polynomial |
---|
256 | za2 = Kb*(1._wp - zba) + K1*(1._wp-zca) |
---|
257 | za1 = K1*Kb*(1._wp - zba - zca) + K1*K2*(1._wp - (zca+zca)) |
---|
258 | za0 = K1*K2*Kb*(1._wp - zba - (zca+zca)) |
---|
259 | ! Taylor expansion around the minimum |
---|
260 | zd = za2*za2 - 3._wp*za1 ! Discriminant of the quadratic equation |
---|
261 | ! for the minimum close to the root |
---|
262 | |
---|
263 | IF(zd > 0._wp) THEN ! If the discriminant is positive |
---|
264 | zsqrtd = SQRT(zd) |
---|
265 | IF(za2 < 0) THEN |
---|
266 | zhmin = (-za2 + zsqrtd)/3._wp |
---|
267 | ELSE |
---|
268 | zhmin = -za1/(za2 + zsqrtd) |
---|
269 | ENDIF |
---|
270 | p_hini = zhmin + SQRT(-(za0 + zhmin*(za1 + zhmin*(za2 + zhmin)))/zsqrtd) |
---|
271 | ELSE |
---|
272 | p_hini = 1.e-7_wp |
---|
273 | ENDIF |
---|
274 | |
---|
275 | ENDIF |
---|
276 | RETURN |
---|
277 | END SUBROUTINE ahini_for_at |
---|
278 | |
---|
279 | !=============================================================================== |
---|
280 | |
---|
281 | FUNCTION solve_at_general(p_alktot, p_dictot, p_bortot, & |
---|
282 | p_po4tot, p_siltot, & |
---|
283 | p_so4tot, p_flutot, & |
---|
284 | K0, K1, K2, Kb, Kw, Ks, Kf, K1p, K2p, K3p, Ksi, & |
---|
285 | p_hini, p_val) |
---|
286 | |
---|
287 | ! Universal pH solver that converges from any given initial value, |
---|
288 | ! determines upper an lower bounds for the solution if required |
---|
289 | |
---|
290 | USE mocsy_singledouble |
---|
291 | IMPLICIT NONE |
---|
292 | REAL(KIND=wp) :: SOLVE_AT_GENERAL |
---|
293 | |
---|
294 | ! Argument variables |
---|
295 | !-------------------- |
---|
296 | REAL(KIND=wp), INTENT(IN) :: p_alktot |
---|
297 | REAL(KIND=wp), INTENT(IN) :: p_dictot |
---|
298 | REAL(KIND=wp), INTENT(IN) :: p_bortot |
---|
299 | REAL(KIND=wp), INTENT(IN) :: p_po4tot |
---|
300 | REAL(KIND=wp), INTENT(IN) :: p_siltot |
---|
301 | !REAL(KIND=wp), INTENT(IN) :: p_nh4tot |
---|
302 | !REAL(KIND=wp), INTENT(IN) :: p_h2stot |
---|
303 | REAL(KIND=wp), INTENT(IN) :: p_so4tot |
---|
304 | REAL(KIND=wp), INTENT(IN) :: p_flutot |
---|
305 | REAL(KIND=wp), INTENT(IN) :: K0, K1, K2, Kb, Kw, Ks, Kf |
---|
306 | REAL(KIND=wp), INTENT(IN) :: K1p, K2p, K3p, Ksi |
---|
307 | REAL(KIND=wp), INTENT(IN), OPTIONAL :: p_hini |
---|
308 | REAL(KIND=wp), INTENT(OUT), OPTIONAL :: p_val |
---|
309 | |
---|
310 | ! Local variables |
---|
311 | !----------------- |
---|
312 | REAL(KIND=wp) :: zh_ini, zh, zh_prev, zh_lnfactor |
---|
313 | REAL(KIND=wp) :: zalknw_inf, zalknw_sup |
---|
314 | REAL(KIND=wp) :: zh_min, zh_max |
---|
315 | REAL(KIND=wp) :: zdelta, zh_delta |
---|
316 | REAL(KIND=wp) :: zeqn, zdeqndh, zeqn_absmin |
---|
317 | REAL(KIND=wp) :: aphscale |
---|
318 | LOGICAL :: l_exitnow |
---|
319 | REAL(KIND=wp), PARAMETER :: pz_exp_threshold = 1.0_wp |
---|
320 | |
---|
321 | ! TOTAL H+ scale: conversion factor for Htot = aphscale * Hfree |
---|
322 | aphscale = 1._wp + p_so4tot/Ks |
---|
323 | |
---|
324 | IF(PRESENT(p_hini)) THEN |
---|
325 | zh_ini = p_hini |
---|
326 | ELSE |
---|
327 | CALL ahini_for_at(p_alktot, p_dictot, p_bortot, K1, K2, Kb, zh_ini) |
---|
328 | ENDIF |
---|
329 | |
---|
330 | CALL anw_infsup(p_dictot, p_bortot, & |
---|
331 | p_po4tot, p_siltot, & |
---|
332 | p_so4tot, p_flutot, & |
---|
333 | zalknw_inf, zalknw_sup) |
---|
334 | |
---|
335 | zdelta = (p_alktot-zalknw_inf)**2 + 4._wp*Kw/aphscale |
---|
336 | |
---|
337 | IF(p_alktot >= zalknw_inf) THEN |
---|
338 | zh_min = 2._wp*Kw /( p_alktot-zalknw_inf + SQRT(zdelta) ) |
---|
339 | ELSE |
---|
340 | zh_min = aphscale*(-(p_alktot-zalknw_inf) + SQRT(zdelta) ) / 2._wp |
---|
341 | ENDIF |
---|
342 | |
---|
343 | zdelta = (p_alktot-zalknw_sup)**2 + 4._wp*Kw/aphscale |
---|
344 | |
---|
345 | IF(p_alktot <= zalknw_sup) THEN |
---|
346 | zh_max = aphscale*(-(p_alktot-zalknw_sup) + SQRT(zdelta) ) / 2._wp |
---|
347 | ELSE |
---|
348 | zh_max = 2._wp*Kw /( p_alktot-zalknw_sup + SQRT(zdelta) ) |
---|
349 | ENDIF |
---|
350 | |
---|
351 | zh = MAX(MIN(zh_max, zh_ini), zh_min) |
---|
352 | niter_atgen = 0 ! Reset counters of iterations |
---|
353 | zeqn_absmin = HUGE(1._wp) |
---|
354 | |
---|
355 | DO |
---|
356 | IF(niter_atgen >= jp_maxniter_atgen) THEN |
---|
357 | zh = -1._wp |
---|
358 | EXIT |
---|
359 | ENDIF |
---|
360 | |
---|
361 | zh_prev = zh |
---|
362 | zeqn = equation_at(p_alktot, zh, p_dictot, p_bortot, & |
---|
363 | p_po4tot, p_siltot, & |
---|
364 | p_so4tot, p_flutot, & |
---|
365 | K0, K1, K2, Kb, Kw, Ks, Kf, K1p, K2p, K3p, Ksi, & |
---|
366 | P_DERIVEQN = zdeqndh) |
---|
367 | |
---|
368 | ! Adapt bracketing interval |
---|
369 | IF(zeqn > 0._wp) THEN |
---|
370 | zh_min = zh_prev |
---|
371 | ELSEIF(zeqn < 0._wp) THEN |
---|
372 | zh_max = zh_prev |
---|
373 | ELSE |
---|
374 | ! zh is the root; unlikely but, one never knows |
---|
375 | EXIT |
---|
376 | ENDIF |
---|
377 | |
---|
378 | ! Now determine the next iterate zh |
---|
379 | niter_atgen = niter_atgen + 1 |
---|
380 | |
---|
381 | IF(ABS(zeqn) >= 0.5_wp*zeqn_absmin) THEN |
---|
382 | ! if the function evaluation at the current point is |
---|
383 | ! not decreasing faster than with a bisection step (at least linearly) |
---|
384 | ! in absolute value take one bisection step on [ph_min, ph_max] |
---|
385 | ! ph_new = (ph_min + ph_max)/2d0 |
---|
386 | ! |
---|
387 | ! In terms of [H]_new: |
---|
388 | ! [H]_new = 10**(-ph_new) |
---|
389 | ! = 10**(-(ph_min + ph_max)/2d0) |
---|
390 | ! = SQRT(10**(-(ph_min + phmax))) |
---|
391 | ! = SQRT(zh_max * zh_min) |
---|
392 | zh = SQRT(zh_max * zh_min) |
---|
393 | zh_lnfactor = (zh - zh_prev)/zh_prev ! Required to test convergence below |
---|
394 | ELSE |
---|
395 | ! dzeqn/dpH = dzeqn/d[H] * d[H]/dpH |
---|
396 | ! = -zdeqndh * LOG(10) * [H] |
---|
397 | ! \Delta pH = -zeqn/(zdeqndh*d[H]/dpH) = zeqn/(zdeqndh*[H]*LOG(10)) |
---|
398 | ! |
---|
399 | ! pH_new = pH_old + \deltapH |
---|
400 | ! |
---|
401 | ! [H]_new = 10**(-pH_new) |
---|
402 | ! = 10**(-pH_old - \Delta pH) |
---|
403 | ! = [H]_old * 10**(-zeqn/(zdeqndh*[H]_old*LOG(10))) |
---|
404 | ! = [H]_old * EXP(-LOG(10)*zeqn/(zdeqndh*[H]_old*LOG(10))) |
---|
405 | ! = [H]_old * EXP(-zeqn/(zdeqndh*[H]_old)) |
---|
406 | |
---|
407 | zh_lnfactor = -zeqn/(zdeqndh*zh_prev) |
---|
408 | |
---|
409 | IF(ABS(zh_lnfactor) > pz_exp_threshold) THEN |
---|
410 | zh = zh_prev*EXP(zh_lnfactor) |
---|
411 | ELSE |
---|
412 | zh_delta = zh_lnfactor*zh_prev |
---|
413 | zh = zh_prev + zh_delta |
---|
414 | ENDIF |
---|
415 | |
---|
416 | IF( zh < zh_min ) THEN |
---|
417 | ! if [H]_new < [H]_min |
---|
418 | ! i.e., if ph_new > ph_max then |
---|
419 | ! take one bisection step on [ph_prev, ph_max] |
---|
420 | ! ph_new = (ph_prev + ph_max)/2d0 |
---|
421 | ! In terms of [H]_new: |
---|
422 | ! [H]_new = 10**(-ph_new) |
---|
423 | ! = 10**(-(ph_prev + ph_max)/2d0) |
---|
424 | ! = SQRT(10**(-(ph_prev + phmax))) |
---|
425 | ! = SQRT([H]_old*10**(-ph_max)) |
---|
426 | ! = SQRT([H]_old * zh_min) |
---|
427 | zh = SQRT(zh_prev * zh_min) |
---|
428 | zh_lnfactor = (zh - zh_prev)/zh_prev ! Required to test convergence below |
---|
429 | ENDIF |
---|
430 | |
---|
431 | IF( zh > zh_max ) THEN |
---|
432 | ! if [H]_new > [H]_max |
---|
433 | ! i.e., if ph_new < ph_min, then |
---|
434 | ! take one bisection step on [ph_min, ph_prev] |
---|
435 | ! ph_new = (ph_prev + ph_min)/2d0 |
---|
436 | ! In terms of [H]_new: |
---|
437 | ! [H]_new = 10**(-ph_new) |
---|
438 | ! = 10**(-(ph_prev + ph_min)/2d0) |
---|
439 | ! = SQRT(10**(-(ph_prev + ph_min))) |
---|
440 | ! = SQRT([H]_old*10**(-ph_min)) |
---|
441 | ! = SQRT([H]_old * zhmax) |
---|
442 | zh = SQRT(zh_prev * zh_max) |
---|
443 | zh_lnfactor = (zh - zh_prev)/zh_prev ! Required to test convergence below |
---|
444 | ENDIF |
---|
445 | ENDIF |
---|
446 | |
---|
447 | zeqn_absmin = MIN( ABS(zeqn), zeqn_absmin) |
---|
448 | |
---|
449 | ! Stop iterations once |\delta{[H]}/[H]| < rdel |
---|
450 | ! <=> |(zh - zh_prev)/zh_prev| = |EXP(-zeqn/(zdeqndh*zh_prev)) -1| < rdel |
---|
451 | ! |EXP(-zeqn/(zdeqndh*zh_prev)) -1| ~ |zeqn/(zdeqndh*zh_prev)| |
---|
452 | |
---|
453 | ! Alternatively: |
---|
454 | ! |\Delta pH| = |zeqn/(zdeqndh*zh_prev*LOG(10))| |
---|
455 | ! ~ 1/LOG(10) * |\Delta [H]|/[H] |
---|
456 | ! < 1/LOG(10) * rdel |
---|
457 | |
---|
458 | ! Hence |zeqn/(zdeqndh*zh)| < rdel |
---|
459 | |
---|
460 | ! rdel <-- pp_rdel_ah_target |
---|
461 | |
---|
462 | l_exitnow = (ABS(zh_lnfactor) < pp_rdel_ah_target) |
---|
463 | |
---|
464 | IF(l_exitnow) EXIT |
---|
465 | ENDDO |
---|
466 | |
---|
467 | solve_at_general = zh |
---|
468 | |
---|
469 | IF(PRESENT(p_val)) THEN |
---|
470 | IF(zh > 0._wp) THEN |
---|
471 | p_val = equation_at(p_alktot, zh, p_dictot, p_bortot, & |
---|
472 | p_po4tot, p_siltot, & |
---|
473 | p_so4tot, p_flutot, & |
---|
474 | K0, K1, K2, Kb, Kw, Ks, Kf, K1p, K2p, K3p, Ksi) |
---|
475 | ELSE |
---|
476 | p_val = HUGE(1._wp) |
---|
477 | ENDIF |
---|
478 | ENDIF |
---|
479 | RETURN |
---|
480 | END FUNCTION solve_at_general |
---|
481 | |
---|
482 | !=============================================================================== |
---|
483 | |
---|
484 | FUNCTION solve_at_general_sec(p_alktot, p_dictot, p_bortot, & |
---|
485 | p_po4tot, p_siltot, & |
---|
486 | p_so4tot, p_flutot, & |
---|
487 | K0, K1, K2, Kb, Kw, Ks, Kf, K1p, K2p, K3p, Ksi, & |
---|
488 | p_hini, p_val) |
---|
489 | |
---|
490 | ! Universal pH solver that converges from any given initial value, |
---|
491 | ! determines upper an lower bounds for the solution if required |
---|
492 | |
---|
493 | !USE MOD_CHEMCONST, ONLY: api1_wat, aphscale |
---|
494 | USE mocsy_singledouble |
---|
495 | IMPLICIT NONE |
---|
496 | REAL(KIND=wp) :: SOLVE_AT_GENERAL_SEC |
---|
497 | |
---|
498 | ! Argument variables |
---|
499 | REAL(KIND=wp), INTENT(IN) :: p_alktot |
---|
500 | REAL(KIND=wp), INTENT(IN) :: p_dictot |
---|
501 | REAL(KIND=wp), INTENT(IN) :: p_bortot |
---|
502 | REAL(KIND=wp), INTENT(IN) :: p_po4tot |
---|
503 | REAL(KIND=wp), INTENT(IN) :: p_siltot |
---|
504 | !REAL(KIND=wp), INTENT(IN) :: p_nh4tot |
---|
505 | !REAL(KIND=wp), INTENT(IN) :: p_h2stot |
---|
506 | REAL(KIND=wp), INTENT(IN) :: p_so4tot |
---|
507 | REAL(KIND=wp), INTENT(IN) :: p_flutot |
---|
508 | REAL(KIND=wp), INTENT(IN) :: K0, K1, K2, Kb, Kw, Ks, Kf |
---|
509 | REAL(KIND=wp), INTENT(IN) :: K1p, K2p, K3p, Ksi |
---|
510 | REAL(KIND=wp), INTENT(IN), OPTIONAL :: p_hini |
---|
511 | REAL(KIND=wp), INTENT(OUT), OPTIONAL :: p_val |
---|
512 | |
---|
513 | ! Local variables |
---|
514 | REAL(KIND=wp) :: zh_ini, zh, zh_1, zh_2, zh_delta |
---|
515 | REAL(KIND=wp) :: zalknw_inf, zalknw_sup |
---|
516 | REAL(KIND=wp) :: zh_min, zh_max |
---|
517 | REAL(KIND=wp) :: zeqn, zeqn_1, zeqn_2, zeqn_absmin |
---|
518 | REAL(KIND=wp) :: zdelta |
---|
519 | REAL(KIND=wp) :: aphscale |
---|
520 | LOGICAL :: l_exitnow |
---|
521 | |
---|
522 | ! TOTAL H+ scale: conversion factor for Htot = aphscale * Hfree |
---|
523 | aphscale = 1._wp + p_so4tot/Ks |
---|
524 | |
---|
525 | IF(PRESENT(p_hini)) THEN |
---|
526 | zh_ini = p_hini |
---|
527 | ELSE |
---|
528 | CALL ahini_for_at(p_alktot, p_dictot, p_bortot, K1, K2, Kb, zh_ini) |
---|
529 | ENDIF |
---|
530 | |
---|
531 | CALL anw_infsup(p_dictot, p_bortot, & |
---|
532 | p_po4tot, p_siltot, & |
---|
533 | p_so4tot, p_flutot, & |
---|
534 | zalknw_inf, zalknw_sup) |
---|
535 | |
---|
536 | zdelta = (p_alktot-zalknw_inf)**2 + 4._wp*Kw/aphscale |
---|
537 | |
---|
538 | IF(p_alktot >= zalknw_inf) THEN |
---|
539 | zh_min = 2._wp*Kw /( p_alktot-zalknw_inf + SQRT(zdelta) ) |
---|
540 | ELSE |
---|
541 | zh_min = aphscale*(-(p_alktot-zalknw_inf) + SQRT(zdelta) ) / 2._wp |
---|
542 | ENDIF |
---|
543 | |
---|
544 | zdelta = (p_alktot-zalknw_sup)**2 + 4._wp*Kw/aphscale |
---|
545 | |
---|
546 | IF(p_alktot <= zalknw_sup) THEN |
---|
547 | zh_max = aphscale*(-(p_alktot-zalknw_sup) + SQRT(zdelta) ) / 2._wp |
---|
548 | ELSE |
---|
549 | zh_max = 2._wp*Kw /( p_alktot-zalknw_sup + SQRT(zdelta) ) |
---|
550 | ENDIF |
---|
551 | |
---|
552 | zh = MAX(MIN(zh_max, zh_ini), zh_min) |
---|
553 | niter_atsec = 0 ! Reset counters of iterations |
---|
554 | |
---|
555 | ! Prepare the secant iterations: two initial (zh, zeqn) pairs are required |
---|
556 | ! We have the starting value, that needs to be completed by the evaluation |
---|
557 | ! of the equation value it produces. |
---|
558 | |
---|
559 | ! Complete the initial value with its equation evaluation |
---|
560 | ! (will take the role of the $n-2$ iterate at the first secant evaluation) |
---|
561 | |
---|
562 | niter_atsec = 0 ! zh_2 is the initial value; |
---|
563 | |
---|
564 | zh_2 = zh |
---|
565 | zeqn_2 = equation_at(p_alktot, zh_2, p_dictot, p_bortot, & |
---|
566 | p_po4tot, p_siltot, & |
---|
567 | p_so4tot, p_flutot, & |
---|
568 | K0, K1, K2, Kb, Kw, Ks, Kf, K1p, K2p, K3p, Ksi) |
---|
569 | |
---|
570 | zeqn_absmin = ABS(zeqn_2) |
---|
571 | |
---|
572 | ! Adapt bracketing interval and heuristically set zh_1 |
---|
573 | IF(zeqn_2 < 0._wp) THEN |
---|
574 | ! If zeqn_2 < 0, then we adjust zh_max: |
---|
575 | ! we can be sure that zh_min < zh_2 < zh_max. |
---|
576 | zh_max = zh_2 |
---|
577 | ! for zh_1, try 25% (0.1 pH units) below the current zh_max, |
---|
578 | ! but stay above SQRT(zh_min*zh_max), which would be equivalent |
---|
579 | ! to a bisection step on [pH@zh_min, pH@zh_max] |
---|
580 | zh_1 = MAX(zh_max/1.25_wp, SQRT(zh_min*zh_max)) |
---|
581 | ELSEIF(zeqn_2 > 0._wp) THEN |
---|
582 | ! If zeqn_2 < 0, then we adjust zh_min: |
---|
583 | ! we can be sure that zh_min < zh_2 < zh_max. |
---|
584 | zh_min = zh_2 |
---|
585 | ! for zh_1, try 25% (0.1 pH units) above the current zh_min, |
---|
586 | ! but stay below SQRT(zh_min*zh_max) which would be equivalent |
---|
587 | ! to a bisection step on [pH@zh_min, pH@zh_max] |
---|
588 | zh_1 = MIN(zh_min*1.25_wp, SQRT(zh_min*zh_max)) |
---|
589 | ELSE ! we have got the root; unlikely, but one never knows |
---|
590 | solve_at_general_sec = zh_2 |
---|
591 | IF(PRESENT(p_val)) p_val = zeqn_2 |
---|
592 | RETURN |
---|
593 | ENDIF |
---|
594 | |
---|
595 | ! We now have the first pair completed (zh_2, zeqn_2). |
---|
596 | ! Define the second one (zh_1, zeqn_1), which is also the first iterate. |
---|
597 | ! zh_1 has already been set above |
---|
598 | niter_atsec = 1 ! Update counter of iterations |
---|
599 | |
---|
600 | zeqn_1 = equation_at(p_alktot, zh_1, p_dictot, p_bortot, & |
---|
601 | p_po4tot, p_siltot, & |
---|
602 | p_so4tot, p_flutot, & |
---|
603 | K0, K1, K2, Kb, Kw, Ks, Kf, K1p, K2p, K3p, Ksi) |
---|
604 | |
---|
605 | ! Adapt bracketing interval: we know that zh_1 <= zh <= zh_max (if zeqn_1 > 0) |
---|
606 | ! or zh_min <= zh <= zh_1 (if zeqn_1 < 0), so this can always be done |
---|
607 | IF(zeqn_1 > 0._wp) THEN |
---|
608 | zh_min = zh_1 |
---|
609 | ELSEIF(zeqn_1 < 0._wp) THEN |
---|
610 | zh_max = zh_1 |
---|
611 | ELSE ! zh_1 is the root |
---|
612 | solve_at_general_sec = zh_1 |
---|
613 | IF(PRESENT(p_val)) p_val = zeqn_1 |
---|
614 | ENDIF |
---|
615 | |
---|
616 | IF(ABS(zeqn_1) > zeqn_absmin) THEN ! Swap zh_2 and zh_1 if ABS(zeqn_2) < ABS(zeqn_1) |
---|
617 | ! so that zh_2 and zh_1 lead to decreasing equation |
---|
618 | ! values (in absolute value) |
---|
619 | zh = zh_1 |
---|
620 | zeqn = zeqn_1 |
---|
621 | zh_1 = zh_2 |
---|
622 | zeqn_1 = zeqn_2 |
---|
623 | zh_2 = zh |
---|
624 | zeqn_2 = zeqn |
---|
625 | ELSE |
---|
626 | zeqn_absmin = ABS(zeqn_1) |
---|
627 | ENDIF |
---|
628 | |
---|
629 | ! Pre-calculate the first secant iterate (this is the second iterate) |
---|
630 | niter_atsec = 2 |
---|
631 | |
---|
632 | zh_delta = -zeqn_1/((zeqn_2-zeqn_1)/(zh_2 - zh_1)) |
---|
633 | zh = zh_1 + zh_delta |
---|
634 | |
---|
635 | ! Make sure that zh_min < zh < zh_max (if not, |
---|
636 | ! bisect around zh_1 which is the best estimate) |
---|
637 | IF (zh > zh_max) THEN ! this can only happen if zh_2 < zh_1 |
---|
638 | ! and zeqn_2 > zeqn_1 > 0 |
---|
639 | zh = SQRT(zh_1*zh_max) |
---|
640 | ENDIF |
---|
641 | |
---|
642 | IF (zh < zh_min) THEN ! this can only happen if zh_2 > zh_1 |
---|
643 | ! and zeqn_2 < zeqn_1 < 0 |
---|
644 | zh = SQRT(zh_1*zh_min) |
---|
645 | ENDIF |
---|
646 | |
---|
647 | DO |
---|
648 | IF(niter_atsec >= jp_maxniter_atsec) THEN |
---|
649 | zh = -1._wp |
---|
650 | EXIT |
---|
651 | ENDIF |
---|
652 | |
---|
653 | zeqn = equation_at(p_alktot, zh, p_dictot, p_bortot, & |
---|
654 | p_po4tot, p_siltot, & |
---|
655 | p_so4tot, p_flutot, & |
---|
656 | K0, K1, K2, Kb, Kw, Ks, Kf, K1p, K2p, K3p, Ksi) |
---|
657 | |
---|
658 | ! Adapt bracketing interval: since initially, zh_min <= zh <= zh_max |
---|
659 | ! we are sure that zh will improve either bracket, depending on the sign |
---|
660 | ! of zeqn |
---|
661 | IF(zeqn > 0._wp) THEN |
---|
662 | zh_min = zh |
---|
663 | ELSEIF(zeqn < 0._wp) THEN |
---|
664 | zh_max = zh |
---|
665 | ELSE |
---|
666 | ! zh is the root |
---|
667 | EXIT |
---|
668 | ENDIF |
---|
669 | |
---|
670 | ! start calculation of next iterate |
---|
671 | niter_atsec = niter_atsec + 1 |
---|
672 | |
---|
673 | zh_2 = zh_1 |
---|
674 | zeqn_2 = zeqn_1 |
---|
675 | zh_1 = zh |
---|
676 | zeqn_1 = zeqn |
---|
677 | |
---|
678 | IF(ABS(zeqn) >= 0.5_wp*zeqn_absmin) THEN |
---|
679 | ! if the function evaluation at the current point |
---|
680 | ! is not decreasing faster in absolute value than |
---|
681 | ! we may expect for a bisection step, then take |
---|
682 | ! one bisection step on [ph_min, ph_max] |
---|
683 | ! ph_new = (ph_min + ph_max)/2d0 |
---|
684 | ! In terms of [H]_new: |
---|
685 | ! [H]_new = 10**(-ph_new) |
---|
686 | ! = 10**(-(ph_min + ph_max)/2d0) |
---|
687 | ! = SQRT(10**(-(ph_min + phmax))) |
---|
688 | ! = SQRT(zh_max * zh_min) |
---|
689 | zh = SQRT(zh_max * zh_min) |
---|
690 | zh_delta = zh - zh_1 |
---|
691 | ELSE |
---|
692 | ! \Delta H = -zeqn_1*(h_2 - h_1)/(zeqn_2 - zeqn_1) |
---|
693 | ! H_new = H_1 + \Delta H |
---|
694 | zh_delta = -zeqn_1/((zeqn_2-zeqn_1)/(zh_2 - zh_1)) |
---|
695 | zh = zh_1 + zh_delta |
---|
696 | |
---|
697 | IF( zh < zh_min ) THEN |
---|
698 | ! if [H]_new < [H]_min |
---|
699 | ! i.e., if ph_new > ph_max then |
---|
700 | ! take one bisection step on [ph_prev, ph_max] |
---|
701 | ! ph_new = (ph_prev + ph_max)/2d0 |
---|
702 | ! In terms of [H]_new: |
---|
703 | ! [H]_new = 10**(-ph_new) |
---|
704 | ! = 10**(-(ph_prev + ph_max)/2d0) |
---|
705 | ! = SQRT(10**(-(ph_prev + phmax))) |
---|
706 | ! = SQRT([H]_old*10**(-ph_max)) |
---|
707 | ! = SQRT([H]_old * zh_min) |
---|
708 | zh = SQRT(zh_1 * zh_min) |
---|
709 | zh_delta = zh - zh_1 |
---|
710 | ENDIF |
---|
711 | |
---|
712 | IF( zh > zh_max ) THEN |
---|
713 | ! if [H]_new > [H]_max |
---|
714 | ! i.e., if ph_new < ph_min, then |
---|
715 | ! take one bisection step on [ph_min, ph_prev] |
---|
716 | ! ph_new = (ph_prev + ph_min)/2d0 |
---|
717 | ! In terms of [H]_new: |
---|
718 | ! [H]_new = 10**(-ph_new) |
---|
719 | ! = 10**(-(ph_prev + ph_min)/2d0) |
---|
720 | ! = SQRT(10**(-(ph_prev + ph_min))) |
---|
721 | ! = SQRT([H]_old*10**(-ph_min)) |
---|
722 | ! = SQRT([H]_old * zhmax) |
---|
723 | zh = SQRT(zh_1 * zh_max) |
---|
724 | zh_delta = zh - zh_1 |
---|
725 | ENDIF |
---|
726 | ENDIF |
---|
727 | |
---|
728 | zeqn_absmin = MIN(ABS(zeqn), zeqn_absmin) |
---|
729 | |
---|
730 | ! Stop iterations once |([H]-[H_1])/[H_1]| < rdel |
---|
731 | l_exitnow = (ABS(zh_delta) < pp_rdel_ah_target*zh_1) |
---|
732 | |
---|
733 | IF(l_exitnow) EXIT |
---|
734 | ENDDO |
---|
735 | |
---|
736 | SOLVE_AT_GENERAL_SEC = zh |
---|
737 | |
---|
738 | IF(PRESENT(p_val)) THEN |
---|
739 | IF(zh > 0._wp) THEN |
---|
740 | p_val = equation_at(p_alktot, zh, p_dictot, p_bortot, & |
---|
741 | p_po4tot, p_siltot, & |
---|
742 | p_so4tot, p_flutot, & |
---|
743 | K0, K1, K2, Kb, Kw, Ks, Kf, K1p, K2p, K3p, Ksi) |
---|
744 | ELSE |
---|
745 | p_val = HUGE(1._wp) |
---|
746 | ENDIF |
---|
747 | ENDIF |
---|
748 | |
---|
749 | RETURN |
---|
750 | END FUNCTION SOLVE_AT_GENERAL_SEC |
---|
751 | |
---|
752 | !=============================================================================== |
---|
753 | |
---|
754 | FUNCTION SOLVE_AT_FAST(p_alktot, p_dictot, p_bortot, & |
---|
755 | p_po4tot, p_siltot, & |
---|
756 | p_so4tot, p_flutot, & |
---|
757 | K0, K1, K2, Kb, Kw, Ks, Kf, K1p, K2p, K3p, Ksi, & |
---|
758 | p_hini, p_val) |
---|
759 | |
---|
760 | ! Fast version of SOLVE_AT_GENERAL, without any bounds checking. |
---|
761 | |
---|
762 | USE mocsy_singledouble |
---|
763 | IMPLICIT NONE |
---|
764 | REAL(KIND=wp) :: SOLVE_AT_FAST |
---|
765 | |
---|
766 | ! Argument variables |
---|
767 | REAL(KIND=wp), INTENT(IN) :: p_alktot |
---|
768 | REAL(KIND=wp), INTENT(IN) :: p_dictot |
---|
769 | REAL(KIND=wp), INTENT(IN) :: p_bortot |
---|
770 | REAL(KIND=wp), INTENT(IN) :: p_po4tot |
---|
771 | REAL(KIND=wp), INTENT(IN) :: p_siltot |
---|
772 | !REAL(KIND=wp), INTENT(IN) :: p_nh4tot |
---|
773 | !REAL(KIND=wp), INTENT(IN) :: p_h2stot |
---|
774 | REAL(KIND=wp), INTENT(IN) :: p_so4tot |
---|
775 | REAL(KIND=wp), INTENT(IN) :: p_flutot |
---|
776 | REAL(KIND=wp), INTENT(IN) :: K0, K1, K2, Kb, Kw, Ks, Kf |
---|
777 | REAL(KIND=wp), INTENT(IN) :: K1p, K2p, K3p, Ksi |
---|
778 | REAL(KIND=wp), INTENT(IN), OPTIONAL :: p_hini |
---|
779 | REAL(KIND=wp), INTENT(OUT), OPTIONAL :: p_val |
---|
780 | |
---|
781 | ! Local variables |
---|
782 | REAL(KIND=wp) :: zh_ini, zh, zh_prev, zh_lnfactor |
---|
783 | REAL(KIND=wp) :: zhdelta |
---|
784 | REAL(KIND=wp) :: zeqn, zdeqndh |
---|
785 | !REAL(KIND=wp) :: aphscale |
---|
786 | LOGICAL :: l_exitnow |
---|
787 | REAL(KIND=wp), PARAMETER :: pz_exp_threshold = 1.0_wp |
---|
788 | |
---|
789 | ! TOTAL H+ scale: conversion factor for Htot = aphscale * Hfree |
---|
790 | !aphscale = 1._wp + p_so4tot/Ks |
---|
791 | |
---|
792 | IF(PRESENT(p_hini)) THEN |
---|
793 | zh_ini = p_hini |
---|
794 | ELSE |
---|
795 | CALL AHINI_FOR_AT(p_alktot, p_dictot, p_bortot, K1, K2, Kb, zh_ini) |
---|
796 | ENDIF |
---|
797 | zh = zh_ini |
---|
798 | |
---|
799 | niter_atfast = 0 ! Reset counters of iterations |
---|
800 | DO |
---|
801 | niter_atfast = niter_atfast + 1 |
---|
802 | IF(niter_atfast > jp_maxniter_atfast) THEN |
---|
803 | zh = -1._wp |
---|
804 | EXIT |
---|
805 | ENDIF |
---|
806 | |
---|
807 | zh_prev = zh |
---|
808 | |
---|
809 | zeqn = equation_at(p_alktot, zh, p_dictot, p_bortot, & |
---|
810 | p_po4tot, p_siltot, & |
---|
811 | p_so4tot, p_flutot, & |
---|
812 | K0, K1, K2, Kb, Kw, Ks, Kf, K1p, K2p, K3p, Ksi, & |
---|
813 | P_DERIVEQN = zdeqndh) |
---|
814 | |
---|
815 | IF(zeqn == 0._wp) EXIT ! zh is the root |
---|
816 | |
---|
817 | zh_lnfactor = -zeqn/(zdeqndh*zh_prev) |
---|
818 | IF(ABS(zh_lnfactor) > pz_exp_threshold) THEN |
---|
819 | zh = zh_prev*EXP(zh_lnfactor) |
---|
820 | ELSE |
---|
821 | zhdelta = zh_lnfactor*zh_prev |
---|
822 | zh = zh_prev + zhdelta |
---|
823 | ENDIF |
---|
824 | |
---|
825 | ! Stop iterations once |\delta{[H]}/[H]| < rdel |
---|
826 | ! <=> |(zh - zh_prev)/zh_prev| = |EXP(-zeqn/(zdeqndh*zh_prev)) -1| < rdel |
---|
827 | ! |EXP(-zeqn/(zdeqndh*zh_prev)) -1| ~ |zeqn/(zdeqndh*zh_prev)| |
---|
828 | |
---|
829 | ! Alternatively: |
---|
830 | ! |\Delta pH| = |zeqn/(zdeqndh*zh_prev*LOG(10))| |
---|
831 | ! ~ 1/LOG(10) * |\Delta [H]|/[H] |
---|
832 | ! < 1/LOG(10) * rdel |
---|
833 | |
---|
834 | ! Hence |zeqn/(zdeqndh*zh)| < rdel |
---|
835 | |
---|
836 | ! rdel <- pp_rdel_ah_target |
---|
837 | |
---|
838 | l_exitnow = (ABS(zh_lnfactor) < pp_rdel_ah_target) |
---|
839 | |
---|
840 | IF(l_exitnow) EXIT |
---|
841 | ENDDO |
---|
842 | |
---|
843 | SOLVE_AT_FAST = zh |
---|
844 | |
---|
845 | IF(PRESENT(p_val)) THEN |
---|
846 | IF(zh > 0._wp) THEN |
---|
847 | p_val = equation_at(p_alktot, zh, p_dictot, p_bortot, & |
---|
848 | p_po4tot, p_siltot, & |
---|
849 | p_so4tot, p_flutot, & |
---|
850 | K0, K1, K2, Kb, Kw, Ks, Kf, K1p, K2p, K3p, Ksi) |
---|
851 | ELSE |
---|
852 | p_val = HUGE(1._wp) |
---|
853 | ENDIF |
---|
854 | ENDIF |
---|
855 | |
---|
856 | RETURN |
---|
857 | END FUNCTION solve_at_fast |
---|
858 | !=============================================================================== |
---|
859 | |
---|
860 | END MODULE mocsy_phsolvers |
---|