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