New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
mocsy_phsolvers.F90 in branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/TOP_SRC/MEDUSA – NEMO

source: branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/TOP_SRC/MEDUSA/mocsy_phsolvers.F90 @ 10020

Last change on this file since 10020 was 5841, checked in by jpalmier, 8 years ago

JPALM --30-10-2015-- Add MOCSY and DMS to MEDUSA-NEMO3.6
File size: 30.4 KB
Line 
1!> \file mocsy_phsolvers.f90
2!! \BRIEF
3!> Module with routines needed to solve pH-total alkalinity equation (Munhoven, 2013, GMD)
4MODULE 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
18USE mocsy_singledouble
19IMPLICIT NONE
20
21! General parameters
22REAL(KIND=wp), PARAMETER :: pp_rdel_ah_target = 1.E-8_wp
23REAL(KIND=wp), PARAMETER :: pp_ln10 = 2.302585092994045684018_wp
24
25! Maximum number of iterations for each method
26INTEGER, PARAMETER :: jp_maxniter_atgen    = 50
27INTEGER, PARAMETER :: jp_maxniter_atsec    = 50
28INTEGER, PARAMETER :: jp_maxniter_atfast   = 50
29
30! Bookkeeping variables for each method
31! - SOLVE_AT_GENERAL
32INTEGER :: niter_atgen    = jp_maxniter_atgen
33
34! - SOLVE_AT_GENERAL_SEC
35INTEGER :: niter_atsec    = jp_maxniter_atsec
36
37! - SOLVE_AT_FAST (variant of SOLVE_AT_GENERAL w/o bracketing
38INTEGER :: 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
44CONTAINS
45!===============================================================================
46SUBROUTINE 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
55USE mocsy_singledouble
56IMPLICIT NONE
57
58! Argument variables
59REAL(KIND=wp), INTENT(IN)  :: p_dictot
60REAL(KIND=wp), INTENT(IN)  :: p_bortot
61REAL(KIND=wp), INTENT(IN)  :: p_po4tot
62REAL(KIND=wp), INTENT(IN)  :: p_siltot
63!REAL(KIND=wp), INTENT(IN)  :: p_nh4tot
64!REAL(KIND=wp), INTENT(IN)  :: p_h2stot
65REAL(KIND=wp), INTENT(IN)  :: p_so4tot
66REAL(KIND=wp), INTENT(IN)  :: p_flutot
67REAL(KIND=wp), INTENT(OUT) :: p_alknw_inf
68REAL(KIND=wp), INTENT(OUT) :: p_alknw_sup
69
70p_alknw_inf =  -p_po4tot - p_so4tot - p_flutot
71p_alknw_sup =   p_dictot + p_dictot + p_bortot &
72              + p_po4tot + p_po4tot + p_siltot !&
73!             + p_nh4tot + p_h2stot
74
75RETURN
76END SUBROUTINE anw_infsup
77
78!===============================================================================
79
80FUNCTION 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
86USE mocsy_singledouble
87IMPLICIT NONE
88REAL(KIND=wp) :: equation_at
89
90! Argument variables
91REAL(KIND=wp), INTENT(IN)            :: p_alktot
92REAL(KIND=wp), INTENT(IN)            :: p_h
93REAL(KIND=wp), INTENT(IN)            :: p_dictot
94REAL(KIND=wp), INTENT(IN)            :: p_bortot
95REAL(KIND=wp), INTENT(IN)            :: p_po4tot
96REAL(KIND=wp), INTENT(IN)            :: p_siltot
97!REAL(KIND=wp), INTENT(IN)            :: p_nh4tot
98!REAL(KIND=wp), INTENT(IN)            :: p_h2stot
99REAL(KIND=wp), INTENT(IN)            :: p_so4tot
100REAL(KIND=wp), INTENT(IN)            :: p_flutot
101REAL(KIND=wp), INTENT(IN)            :: K0, K1, K2, Kb, Kw, Ks, Kf
102REAL(KIND=wp), INTENT(IN)            :: K1p, K2p, K3p, Ksi
103REAL(KIND=wp), INTENT(OUT), OPTIONAL :: p_deriveqn
104
105! Local variables
106!-----------------
107REAL(KIND=wp) :: znumer_dic, zdnumer_dic, zdenom_dic, zalk_dic, zdalk_dic
108REAL(KIND=wp) :: znumer_bor, zdnumer_bor, zdenom_bor, zalk_bor, zdalk_bor
109REAL(KIND=wp) :: znumer_po4, zdnumer_po4, zdenom_po4, zalk_po4, zdalk_po4
110REAL(KIND=wp) :: znumer_sil, zdnumer_sil, zdenom_sil, zalk_sil, zdalk_sil
111REAL(KIND=wp) :: znumer_nh4, zdnumer_nh4, zdenom_nh4, zalk_nh4, zdalk_nh4
112REAL(KIND=wp) :: znumer_h2s, zdnumer_h2s, zdenom_h2s, zalk_h2s, zdalk_h2s
113REAL(KIND=wp) :: znumer_so4, zdnumer_so4, zdenom_so4, zalk_so4, zdalk_so4
114REAL(KIND=wp) :: znumer_flu, zdnumer_flu, zdenom_flu, zalk_flu, zdalk_flu
115REAL(KIND=wp) ::                                      zalk_wat, zdalk_wat
116REAL(KIND=wp) :: aphscale
117
118! TOTAL H+ scale: conversion factor for Htot = aphscale * Hfree
119aphscale = 1._wp + p_so4tot/Ks
120
121! H2CO3 - HCO3 - CO3 : n=2, m=0
122znumer_dic = 2._wp*K1*K2 + p_h*       K1
123zdenom_dic =       K1*K2 + p_h*(      K1 + p_h)
124zalk_dic   = p_dictot * (znumer_dic/zdenom_dic)
125
126! B(OH)3 - B(OH)4 : n=1, m=0
127znumer_bor =       Kb
128zdenom_bor =       Kb + p_h
129zalk_bor   = p_bortot * (znumer_bor/zdenom_bor)
130
131! H3PO4 - H2PO4 - HPO4 - PO4 : n=3, m=1
132znumer_po4 = 3._wp*K1p*K2p*K3p + p_h*(2._wp*K1p*K2p + p_h* K1p)
133zdenom_po4 =       K1p*K2p*K3p + p_h*(      K1p*K2p + p_h*(K1p + p_h))
134zalk_po4   = p_po4tot * (znumer_po4/zdenom_po4 - 1._wp) ! Zero level of H3PO4 = 1
135
136! H4SiO4 - H3SiO4 : n=1, m=0
137znumer_sil =       Ksi
138zdenom_sil =       Ksi + p_h
139zalk_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
154znumer_so4 =       Ks
155zdenom_so4 =       Ks + p_h
156zalk_so4   = p_so4tot * (znumer_so4/zdenom_so4 - 1._wp)
157
158! HF - F : n=1, m=1
159znumer_flu =       Kf
160zdenom_flu =       Kf + p_h
161zalk_flu   = p_flutot * (znumer_flu/zdenom_flu - 1._wp)
162
163! H2O - OH
164zalk_wat   = Kw/p_h - p_h/aphscale
165
166equation_at =    zalk_dic + zalk_bor + zalk_po4 + zalk_sil &
167               + zalk_so4 + zalk_flu                       &
168               + zalk_wat - p_alktot
169
170IF(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
213ENDIF
214RETURN
215END FUNCTION equation_at
216
217!===============================================================================
218
219SUBROUTINE 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
232USE mocsy_singledouble
233IMPLICIT NONE
234
235! Argument variables
236!--------------------
237REAL(KIND=wp), INTENT(IN)   ::  p_alkcb, p_dictot, p_bortot
238REAL(KIND=wp), INTENT(IN)   ::  K1, K2, Kb
239REAL(KIND=wp), INTENT(OUT)  ::  p_hini
240
241! Local variables
242!-----------------
243REAL(KIND=wp)  ::  zca, zba
244REAL(KIND=wp)  ::  zd, zsqrtd, zhmin
245REAL(KIND=wp)  ::  za2, za1, za0
246
247IF (p_alkcb <= 0._wp) THEN
248  p_hini = 1.e-3_wp
249ELSEIF (p_alkcb >= (2._wp*p_dictot + p_bortot)) THEN
250  p_hini = 1.e-10_wp
251ELSE
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
275ENDIF
276RETURN
277END SUBROUTINE ahini_for_at
278
279!===============================================================================
280
281FUNCTION 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
290USE mocsy_singledouble
291IMPLICIT NONE
292REAL(KIND=wp) :: SOLVE_AT_GENERAL
293
294! Argument variables
295!--------------------
296REAL(KIND=wp), INTENT(IN)            :: p_alktot
297REAL(KIND=wp), INTENT(IN)            :: p_dictot
298REAL(KIND=wp), INTENT(IN)            :: p_bortot
299REAL(KIND=wp), INTENT(IN)            :: p_po4tot
300REAL(KIND=wp), INTENT(IN)            :: p_siltot
301!REAL(KIND=wp), INTENT(IN)            :: p_nh4tot
302!REAL(KIND=wp), INTENT(IN)            :: p_h2stot
303REAL(KIND=wp), INTENT(IN)            :: p_so4tot
304REAL(KIND=wp), INTENT(IN)            :: p_flutot
305REAL(KIND=wp), INTENT(IN)            :: K0, K1, K2, Kb, Kw, Ks, Kf
306REAL(KIND=wp), INTENT(IN)            :: K1p, K2p, K3p, Ksi
307REAL(KIND=wp), INTENT(IN), OPTIONAL  :: p_hini
308REAL(KIND=wp), INTENT(OUT), OPTIONAL :: p_val
309
310! Local variables
311!-----------------
312REAL(KIND=wp)  ::  zh_ini, zh, zh_prev, zh_lnfactor
313REAL(KIND=wp)  ::  zalknw_inf, zalknw_sup
314REAL(KIND=wp)  ::  zh_min, zh_max
315REAL(KIND=wp)  ::  zdelta, zh_delta
316REAL(KIND=wp)  ::  zeqn, zdeqndh, zeqn_absmin
317REAL(KIND=wp)  ::  aphscale
318LOGICAL        :: l_exitnow
319REAL(KIND=wp), PARAMETER :: pz_exp_threshold = 1.0_wp
320
321! TOTAL H+ scale: conversion factor for Htot = aphscale * Hfree
322aphscale = 1._wp + p_so4tot/Ks
323
324IF(PRESENT(p_hini)) THEN
325   zh_ini = p_hini
326ELSE
327   CALL ahini_for_at(p_alktot, p_dictot, p_bortot, K1, K2, Kb, zh_ini)
328ENDIF
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
335zdelta = (p_alktot-zalknw_inf)**2 + 4._wp*Kw/aphscale
336
337IF(p_alktot >= zalknw_inf) THEN
338   zh_min = 2._wp*Kw /( p_alktot-zalknw_inf + SQRT(zdelta) )
339ELSE
340   zh_min = aphscale*(-(p_alktot-zalknw_inf) + SQRT(zdelta) ) / 2._wp
341ENDIF
342
343zdelta = (p_alktot-zalknw_sup)**2 + 4._wp*Kw/aphscale
344
345IF(p_alktot <= zalknw_sup) THEN
346   zh_max = aphscale*(-(p_alktot-zalknw_sup) + SQRT(zdelta) ) / 2._wp
347ELSE
348   zh_max = 2._wp*Kw /( p_alktot-zalknw_sup + SQRT(zdelta) )
349ENDIF
350
351zh = MAX(MIN(zh_max, zh_ini), zh_min)
352niter_atgen        = 0                 ! Reset counters of iterations
353zeqn_absmin        = HUGE(1._wp)
354
355DO
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
465ENDDO
466
467solve_at_general = zh
468
469IF(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
478ENDIF
479RETURN
480END FUNCTION solve_at_general
481
482!===============================================================================
483
484FUNCTION 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
494USE mocsy_singledouble
495IMPLICIT NONE
496REAL(KIND=wp) :: SOLVE_AT_GENERAL_SEC
497
498! Argument variables
499REAL(KIND=wp), INTENT(IN)            :: p_alktot
500REAL(KIND=wp), INTENT(IN)            :: p_dictot
501REAL(KIND=wp), INTENT(IN)            :: p_bortot
502REAL(KIND=wp), INTENT(IN)            :: p_po4tot
503REAL(KIND=wp), INTENT(IN)            :: p_siltot
504!REAL(KIND=wp), INTENT(IN)            :: p_nh4tot
505!REAL(KIND=wp), INTENT(IN)            :: p_h2stot
506REAL(KIND=wp), INTENT(IN)            :: p_so4tot
507REAL(KIND=wp), INTENT(IN)            :: p_flutot
508REAL(KIND=wp), INTENT(IN)            :: K0, K1, K2, Kb, Kw, Ks, Kf
509REAL(KIND=wp), INTENT(IN)            :: K1p, K2p, K3p, Ksi
510REAL(KIND=wp), INTENT(IN), OPTIONAL  :: p_hini
511REAL(KIND=wp), INTENT(OUT), OPTIONAL :: p_val
512
513! Local variables
514REAL(KIND=wp)  ::  zh_ini, zh, zh_1, zh_2, zh_delta
515REAL(KIND=wp)  ::  zalknw_inf, zalknw_sup
516REAL(KIND=wp)  ::  zh_min, zh_max
517REAL(KIND=wp)  ::  zeqn, zeqn_1, zeqn_2, zeqn_absmin
518REAL(KIND=wp)  ::  zdelta
519REAL(KIND=wp)  ::  aphscale
520LOGICAL        ::  l_exitnow
521
522! TOTAL H+ scale: conversion factor for Htot = aphscale * Hfree
523aphscale = 1._wp + p_so4tot/Ks
524
525IF(PRESENT(p_hini)) THEN
526   zh_ini = p_hini
527ELSE
528   CALL ahini_for_at(p_alktot, p_dictot, p_bortot, K1, K2, Kb, zh_ini)
529ENDIF
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
536zdelta = (p_alktot-zalknw_inf)**2 + 4._wp*Kw/aphscale
537
538IF(p_alktot >= zalknw_inf) THEN
539   zh_min = 2._wp*Kw /( p_alktot-zalknw_inf + SQRT(zdelta) )
540ELSE
541   zh_min = aphscale*(-(p_alktot-zalknw_inf) + SQRT(zdelta) ) / 2._wp
542ENDIF
543
544zdelta = (p_alktot-zalknw_sup)**2 + 4._wp*Kw/aphscale
545
546IF(p_alktot <= zalknw_sup) THEN
547   zh_max = aphscale*(-(p_alktot-zalknw_sup) + SQRT(zdelta) ) / 2._wp
548ELSE
549   zh_max = 2._wp*Kw /( p_alktot-zalknw_sup + SQRT(zdelta) )
550ENDIF
551
552zh = MAX(MIN(zh_max, zh_ini), zh_min)
553niter_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
562niter_atsec = 0                        ! zh_2 is the initial value;
563
564zh_2   = zh
565zeqn_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
570zeqn_absmin        = ABS(zeqn_2)
571
572! Adapt bracketing interval and heuristically set zh_1
573IF(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))
581ELSEIF(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))
589ELSE ! 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
593ENDIF
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
598niter_atsec = 1                        ! Update counter of iterations
599
600zeqn_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
607IF(zeqn_1 > 0._wp) THEN
608   zh_min = zh_1
609ELSEIF(zeqn_1 < 0._wp) THEN
610   zh_max = zh_1
611ELSE ! zh_1 is the root
612   solve_at_general_sec = zh_1
613   IF(PRESENT(p_val)) p_val = zeqn_1
614ENDIF
615
616IF(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
625ELSE
626   zeqn_absmin = ABS(zeqn_1)
627ENDIF
628
629! Pre-calculate the first secant iterate (this is the second iterate)
630niter_atsec = 2
631
632zh_delta = -zeqn_1/((zeqn_2-zeqn_1)/(zh_2 - zh_1))
633zh = 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)
637IF (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)
640ENDIF
641
642IF (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)
645ENDIF
646
647DO
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
734ENDDO
735
736SOLVE_AT_GENERAL_SEC = zh
737
738IF(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
747ENDIF
748
749RETURN
750END FUNCTION SOLVE_AT_GENERAL_SEC
751
752!===============================================================================
753
754FUNCTION 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
762USE mocsy_singledouble
763IMPLICIT NONE
764REAL(KIND=wp) :: SOLVE_AT_FAST
765
766! Argument variables
767REAL(KIND=wp), INTENT(IN)            :: p_alktot
768REAL(KIND=wp), INTENT(IN)            :: p_dictot
769REAL(KIND=wp), INTENT(IN)            :: p_bortot
770REAL(KIND=wp), INTENT(IN)            :: p_po4tot
771REAL(KIND=wp), INTENT(IN)            :: p_siltot
772!REAL(KIND=wp), INTENT(IN)            :: p_nh4tot
773!REAL(KIND=wp), INTENT(IN)            :: p_h2stot
774REAL(KIND=wp), INTENT(IN)            :: p_so4tot
775REAL(KIND=wp), INTENT(IN)            :: p_flutot
776REAL(KIND=wp), INTENT(IN)            :: K0, K1, K2, Kb, Kw, Ks, Kf
777REAL(KIND=wp), INTENT(IN)            :: K1p, K2p, K3p, Ksi
778REAL(KIND=wp), INTENT(IN), OPTIONAL  :: p_hini
779REAL(KIND=wp), INTENT(OUT), OPTIONAL :: p_val
780
781! Local variables
782REAL(KIND=wp)  ::  zh_ini, zh, zh_prev, zh_lnfactor
783REAL(KIND=wp)  ::  zhdelta
784REAL(KIND=wp)  ::  zeqn, zdeqndh
785!REAL(KIND=wp)  ::  aphscale
786LOGICAL        :: l_exitnow
787REAL(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
792IF(PRESENT(p_hini)) THEN
793   zh_ini = p_hini
794ELSE
795   CALL AHINI_FOR_AT(p_alktot, p_dictot, p_bortot, K1, K2, Kb, zh_ini)
796ENDIF
797zh = zh_ini
798
799niter_atfast    = 0                 ! Reset counters of iterations
800DO
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
841ENDDO
842
843SOLVE_AT_FAST = zh
844
845IF(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
854ENDIF
855
856RETURN
857END FUNCTION solve_at_fast
858!===============================================================================
859
860END MODULE mocsy_phsolvers
Note: See TracBrowser for help on using the repository browser.