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_under_ice_relax_dr_hook/NEMOGCM/NEMO/TOP_SRC/MEDUSA – NEMO

source: branches/UKMO/dev_r5518_GO6_under_ice_relax_dr_hook/NEMOGCM/NEMO/TOP_SRC/MEDUSA/mocsy_phsolvers.F90 @ 11738

Last change on this file since 11738 was 11738, checked in by marc, 4 years ago

The Dr Hook changes from my perl code.

File size: 33.2 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
19USE yomhook, ONLY: lhook, dr_hook
20USE parkind1, ONLY: jprb, jpim
21
22IMPLICIT NONE
23
24! General parameters
25REAL(KIND=wp), PARAMETER :: pp_rdel_ah_target = 1.E-8_wp
26REAL(KIND=wp), PARAMETER :: pp_ln10 = 2.302585092994045684018_wp
27
28! Maximum number of iterations for each method
29INTEGER, PARAMETER :: jp_maxniter_atgen    = 50
30INTEGER, PARAMETER :: jp_maxniter_atsec    = 50
31INTEGER, PARAMETER :: jp_maxniter_atfast   = 50
32
33! Bookkeeping variables for each method
34! - SOLVE_AT_GENERAL
35INTEGER :: niter_atgen    = jp_maxniter_atgen
36
37! - SOLVE_AT_GENERAL_SEC
38INTEGER :: niter_atsec    = jp_maxniter_atsec
39
40! - SOLVE_AT_FAST (variant of SOLVE_AT_GENERAL w/o bracketing
41INTEGER :: 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
47CONTAINS
48!===============================================================================
49SUBROUTINE 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
58USE mocsy_singledouble
59USE yomhook, ONLY: lhook, dr_hook
60USE parkind1, ONLY: jprb, jpim
61
62IMPLICIT NONE
63
64! Argument variables
65REAL(KIND=wp), INTENT(IN)  :: p_dictot
66REAL(KIND=wp), INTENT(IN)  :: p_bortot
67REAL(KIND=wp), INTENT(IN)  :: p_po4tot
68REAL(KIND=wp), INTENT(IN)  :: p_siltot
69!REAL(KIND=wp), INTENT(IN)  :: p_nh4tot
70!REAL(KIND=wp), INTENT(IN)  :: p_h2stot
71REAL(KIND=wp), INTENT(IN)  :: p_so4tot
72REAL(KIND=wp), INTENT(IN)  :: p_flutot
73REAL(KIND=wp), INTENT(OUT) :: p_alknw_inf
74REAL(KIND=wp), INTENT(OUT) :: p_alknw_sup
75INTEGER(KIND=jpim), PARAMETER :: zhook_in = 0
76INTEGER(KIND=jpim), PARAMETER :: zhook_out = 1
77REAL(KIND=jprb)               :: zhook_handle
78
79CHARACTER(LEN=*), PARAMETER :: RoutineName='ANW_INFSUP'
80
81IF (lhook) CALL dr_hook(RoutineName,zhook_in,zhook_handle)
82
83
84p_alknw_inf =  -p_po4tot - p_so4tot - p_flutot
85p_alknw_sup =   p_dictot + p_dictot + p_bortot &
86              + p_po4tot + p_po4tot + p_siltot !&
87!             + p_nh4tot + p_h2stot
88
89IF (lhook) CALL dr_hook(RoutineName,zhook_out ,zhook_handle)
90RETURN
91IF (lhook) CALL dr_hook(RoutineName,zhook_out,zhook_handle)
92END SUBROUTINE anw_infsup
93
94!===============================================================================
95
96FUNCTION 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
102USE mocsy_singledouble
103USE yomhook, ONLY: lhook, dr_hook
104USE parkind1, ONLY: jprb, jpim
105
106IMPLICIT NONE
107REAL(KIND=wp) :: equation_at
108
109! Argument variables
110REAL(KIND=wp), INTENT(IN)            :: p_alktot
111REAL(KIND=wp), INTENT(IN)            :: p_h
112REAL(KIND=wp), INTENT(IN)            :: p_dictot
113REAL(KIND=wp), INTENT(IN)            :: p_bortot
114REAL(KIND=wp), INTENT(IN)            :: p_po4tot
115REAL(KIND=wp), INTENT(IN)            :: p_siltot
116!REAL(KIND=wp), INTENT(IN)            :: p_nh4tot
117!REAL(KIND=wp), INTENT(IN)            :: p_h2stot
118REAL(KIND=wp), INTENT(IN)            :: p_so4tot
119REAL(KIND=wp), INTENT(IN)            :: p_flutot
120REAL(KIND=wp), INTENT(IN)            :: K0, K1, K2, Kb, Kw, Ks, Kf
121REAL(KIND=wp), INTENT(IN)            :: K1p, K2p, K3p, Ksi
122REAL(KIND=wp), INTENT(OUT), OPTIONAL :: p_deriveqn
123
124! Local variables
125!-----------------
126REAL(KIND=wp) :: znumer_dic, zdnumer_dic, zdenom_dic, zalk_dic, zdalk_dic
127REAL(KIND=wp) :: znumer_bor, zdnumer_bor, zdenom_bor, zalk_bor, zdalk_bor
128REAL(KIND=wp) :: znumer_po4, zdnumer_po4, zdenom_po4, zalk_po4, zdalk_po4
129REAL(KIND=wp) :: znumer_sil, zdnumer_sil, zdenom_sil, zalk_sil, zdalk_sil
130REAL(KIND=wp) :: znumer_nh4, zdnumer_nh4, zdenom_nh4, zalk_nh4, zdalk_nh4
131REAL(KIND=wp) :: znumer_h2s, zdnumer_h2s, zdenom_h2s, zalk_h2s, zdalk_h2s
132REAL(KIND=wp) :: znumer_so4, zdnumer_so4, zdenom_so4, zalk_so4, zdalk_so4
133REAL(KIND=wp) :: znumer_flu, zdnumer_flu, zdenom_flu, zalk_flu, zdalk_flu
134REAL(KIND=wp) ::                                      zalk_wat, zdalk_wat
135REAL(KIND=wp) :: aphscale
136INTEGER(KIND=jpim), PARAMETER :: zhook_in = 0
137INTEGER(KIND=jpim), PARAMETER :: zhook_out = 1
138REAL(KIND=jprb)               :: zhook_handle
139
140CHARACTER(LEN=*), PARAMETER :: RoutineName='EQUATION_AT'
141
142IF (lhook) CALL dr_hook(RoutineName,zhook_in,zhook_handle)
143
144
145! TOTAL H+ scale: conversion factor for Htot = aphscale * Hfree
146aphscale = 1._wp + p_so4tot/Ks
147
148! H2CO3 - HCO3 - CO3 : n=2, m=0
149znumer_dic = 2._wp*K1*K2 + p_h*       K1
150zdenom_dic =       K1*K2 + p_h*(      K1 + p_h)
151zalk_dic   = p_dictot * (znumer_dic/zdenom_dic)
152
153! B(OH)3 - B(OH)4 : n=1, m=0
154znumer_bor =       Kb
155zdenom_bor =       Kb + p_h
156zalk_bor   = p_bortot * (znumer_bor/zdenom_bor)
157
158! H3PO4 - H2PO4 - HPO4 - PO4 : n=3, m=1
159znumer_po4 = 3._wp*K1p*K2p*K3p + p_h*(2._wp*K1p*K2p + p_h* K1p)
160zdenom_po4 =       K1p*K2p*K3p + p_h*(      K1p*K2p + p_h*(K1p + p_h))
161zalk_po4   = p_po4tot * (znumer_po4/zdenom_po4 - 1._wp) ! Zero level of H3PO4 = 1
162
163! H4SiO4 - H3SiO4 : n=1, m=0
164znumer_sil =       Ksi
165zdenom_sil =       Ksi + p_h
166zalk_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
181znumer_so4 =       Ks
182zdenom_so4 =       Ks + p_h
183zalk_so4   = p_so4tot * (znumer_so4/zdenom_so4 - 1._wp)
184
185! HF - F : n=1, m=1
186znumer_flu =       Kf
187zdenom_flu =       Kf + p_h
188zalk_flu   = p_flutot * (znumer_flu/zdenom_flu - 1._wp)
189
190! H2O - OH
191zalk_wat   = Kw/p_h - p_h/aphscale
192
193equation_at =    zalk_dic + zalk_bor + zalk_po4 + zalk_sil &
194               + zalk_so4 + zalk_flu                       &
195               + zalk_wat - p_alktot
196
197IF(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
240ENDIF
241IF (lhook) CALL dr_hook(RoutineName,zhook_out ,zhook_handle)
242RETURN
243IF (lhook) CALL dr_hook(RoutineName,zhook_out,zhook_handle)
244END FUNCTION equation_at
245
246!===============================================================================
247
248SUBROUTINE 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
261USE mocsy_singledouble
262USE yomhook, ONLY: lhook, dr_hook
263USE parkind1, ONLY: jprb, jpim
264
265IMPLICIT NONE
266
267! Argument variables
268!--------------------
269REAL(KIND=wp), INTENT(IN)   ::  p_alkcb, p_dictot, p_bortot
270REAL(KIND=wp), INTENT(IN)   ::  K1, K2, Kb
271REAL(KIND=wp), INTENT(OUT)  ::  p_hini
272
273! Local variables
274!-----------------
275REAL(KIND=wp)  ::  zca, zba
276REAL(KIND=wp)  ::  zd, zsqrtd, zhmin
277REAL(KIND=wp)  ::  za2, za1, za0
278INTEGER(KIND=jpim), PARAMETER :: zhook_in = 0
279INTEGER(KIND=jpim), PARAMETER :: zhook_out = 1
280REAL(KIND=jprb)               :: zhook_handle
281
282CHARACTER(LEN=*), PARAMETER :: RoutineName='AHINI_FOR_AT'
283
284IF (lhook) CALL dr_hook(RoutineName,zhook_in,zhook_handle)
285
286
287IF (p_alkcb <= 0._wp) THEN
288  p_hini = 1.e-3_wp
289ELSEIF (p_alkcb >= (2._wp*p_dictot + p_bortot)) THEN
290  p_hini = 1.e-10_wp
291ELSE
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
315ENDIF
316IF (lhook) CALL dr_hook(RoutineName,zhook_out ,zhook_handle)
317RETURN
318IF (lhook) CALL dr_hook(RoutineName,zhook_out,zhook_handle)
319END SUBROUTINE ahini_for_at
320
321!===============================================================================
322
323FUNCTION 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
332USE mocsy_singledouble
333USE yomhook, ONLY: lhook, dr_hook
334USE parkind1, ONLY: jprb, jpim
335
336IMPLICIT NONE
337REAL(KIND=wp) :: SOLVE_AT_GENERAL
338
339! Argument variables
340!--------------------
341REAL(KIND=wp), INTENT(IN)            :: p_alktot
342REAL(KIND=wp), INTENT(IN)            :: p_dictot
343REAL(KIND=wp), INTENT(IN)            :: p_bortot
344REAL(KIND=wp), INTENT(IN)            :: p_po4tot
345REAL(KIND=wp), INTENT(IN)            :: p_siltot
346!REAL(KIND=wp), INTENT(IN)            :: p_nh4tot
347!REAL(KIND=wp), INTENT(IN)            :: p_h2stot
348REAL(KIND=wp), INTENT(IN)            :: p_so4tot
349REAL(KIND=wp), INTENT(IN)            :: p_flutot
350REAL(KIND=wp), INTENT(IN)            :: K0, K1, K2, Kb, Kw, Ks, Kf
351REAL(KIND=wp), INTENT(IN)            :: K1p, K2p, K3p, Ksi
352REAL(KIND=wp), INTENT(IN), OPTIONAL  :: p_hini
353REAL(KIND=wp), INTENT(OUT), OPTIONAL :: p_val
354
355! Local variables
356!-----------------
357REAL(KIND=wp)  ::  zh_ini, zh, zh_prev, zh_lnfactor
358REAL(KIND=wp)  ::  zalknw_inf, zalknw_sup
359REAL(KIND=wp)  ::  zh_min, zh_max
360REAL(KIND=wp)  ::  zdelta, zh_delta
361REAL(KIND=wp)  ::  zeqn, zdeqndh, zeqn_absmin
362REAL(KIND=wp)  ::  aphscale
363LOGICAL        :: l_exitnow
364REAL(KIND=wp), PARAMETER :: pz_exp_threshold = 1.0_wp
365INTEGER(KIND=jpim), PARAMETER :: zhook_in = 0
366INTEGER(KIND=jpim), PARAMETER :: zhook_out = 1
367REAL(KIND=jprb)               :: zhook_handle
368
369CHARACTER(LEN=*), PARAMETER :: RoutineName='SOLVE_AT_GENERAL'
370
371IF (lhook) CALL dr_hook(RoutineName,zhook_in,zhook_handle)
372
373
374! TOTAL H+ scale: conversion factor for Htot = aphscale * Hfree
375aphscale = 1._wp + p_so4tot/Ks
376
377IF(PRESENT(p_hini)) THEN
378   zh_ini = p_hini
379ELSE
380   CALL ahini_for_at(p_alktot, p_dictot, p_bortot, K1, K2, Kb, zh_ini)
381ENDIF
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
388zdelta = (p_alktot-zalknw_inf)**2 + 4._wp*Kw/aphscale
389
390IF(p_alktot >= zalknw_inf) THEN
391   zh_min = 2._wp*Kw /( p_alktot-zalknw_inf + SQRT(zdelta) )
392ELSE
393   zh_min = aphscale*(-(p_alktot-zalknw_inf) + SQRT(zdelta) ) / 2._wp
394ENDIF
395
396zdelta = (p_alktot-zalknw_sup)**2 + 4._wp*Kw/aphscale
397
398IF(p_alktot <= zalknw_sup) THEN
399   zh_max = aphscale*(-(p_alktot-zalknw_sup) + SQRT(zdelta) ) / 2._wp
400ELSE
401   zh_max = 2._wp*Kw /( p_alktot-zalknw_sup + SQRT(zdelta) )
402ENDIF
403
404zh = MAX(MIN(zh_max, zh_ini), zh_min)
405niter_atgen        = 0                 ! Reset counters of iterations
406zeqn_absmin        = HUGE(1._wp)
407
408DO
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
518ENDDO
519
520solve_at_general = zh
521
522IF(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
531ENDIF
532IF (lhook) CALL dr_hook(RoutineName,zhook_out ,zhook_handle)
533RETURN
534IF (lhook) CALL dr_hook(RoutineName,zhook_out,zhook_handle)
535END FUNCTION solve_at_general
536
537!===============================================================================
538
539FUNCTION 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
549USE mocsy_singledouble
550USE yomhook, ONLY: lhook, dr_hook
551USE parkind1, ONLY: jprb, jpim
552
553IMPLICIT NONE
554REAL(KIND=wp) :: SOLVE_AT_GENERAL_SEC
555
556! Argument variables
557REAL(KIND=wp), INTENT(IN)            :: p_alktot
558REAL(KIND=wp), INTENT(IN)            :: p_dictot
559REAL(KIND=wp), INTENT(IN)            :: p_bortot
560REAL(KIND=wp), INTENT(IN)            :: p_po4tot
561REAL(KIND=wp), INTENT(IN)            :: p_siltot
562!REAL(KIND=wp), INTENT(IN)            :: p_nh4tot
563!REAL(KIND=wp), INTENT(IN)            :: p_h2stot
564REAL(KIND=wp), INTENT(IN)            :: p_so4tot
565REAL(KIND=wp), INTENT(IN)            :: p_flutot
566REAL(KIND=wp), INTENT(IN)            :: K0, K1, K2, Kb, Kw, Ks, Kf
567REAL(KIND=wp), INTENT(IN)            :: K1p, K2p, K3p, Ksi
568REAL(KIND=wp), INTENT(IN), OPTIONAL  :: p_hini
569REAL(KIND=wp), INTENT(OUT), OPTIONAL :: p_val
570
571! Local variables
572REAL(KIND=wp)  ::  zh_ini, zh, zh_1, zh_2, zh_delta
573REAL(KIND=wp)  ::  zalknw_inf, zalknw_sup
574REAL(KIND=wp)  ::  zh_min, zh_max
575REAL(KIND=wp)  ::  zeqn, zeqn_1, zeqn_2, zeqn_absmin
576REAL(KIND=wp)  ::  zdelta
577REAL(KIND=wp)  ::  aphscale
578LOGICAL        ::  l_exitnow
579INTEGER(KIND=jpim), PARAMETER :: zhook_in = 0
580INTEGER(KIND=jpim), PARAMETER :: zhook_out = 1
581REAL(KIND=jprb)               :: zhook_handle
582
583CHARACTER(LEN=*), PARAMETER :: RoutineName='SOLVE_AT_GENERAL_SEC'
584
585IF (lhook) CALL dr_hook(RoutineName,zhook_in,zhook_handle)
586
587
588! TOTAL H+ scale: conversion factor for Htot = aphscale * Hfree
589aphscale = 1._wp + p_so4tot/Ks
590
591IF(PRESENT(p_hini)) THEN
592   zh_ini = p_hini
593ELSE
594   CALL ahini_for_at(p_alktot, p_dictot, p_bortot, K1, K2, Kb, zh_ini)
595ENDIF
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
602zdelta = (p_alktot-zalknw_inf)**2 + 4._wp*Kw/aphscale
603
604IF(p_alktot >= zalknw_inf) THEN
605   zh_min = 2._wp*Kw /( p_alktot-zalknw_inf + SQRT(zdelta) )
606ELSE
607   zh_min = aphscale*(-(p_alktot-zalknw_inf) + SQRT(zdelta) ) / 2._wp
608ENDIF
609
610zdelta = (p_alktot-zalknw_sup)**2 + 4._wp*Kw/aphscale
611
612IF(p_alktot <= zalknw_sup) THEN
613   zh_max = aphscale*(-(p_alktot-zalknw_sup) + SQRT(zdelta) ) / 2._wp
614ELSE
615   zh_max = 2._wp*Kw /( p_alktot-zalknw_sup + SQRT(zdelta) )
616ENDIF
617
618zh = MAX(MIN(zh_max, zh_ini), zh_min)
619niter_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
628niter_atsec = 0                        ! zh_2 is the initial value;
629
630zh_2   = zh
631zeqn_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
636zeqn_absmin        = ABS(zeqn_2)
637
638! Adapt bracketing interval and heuristically set zh_1
639IF(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))
647ELSEIF(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))
655ELSE ! 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
658IF (lhook) CALL dr_hook(RoutineName,zhook_out ,zhook_handle)
659   RETURN
660ENDIF
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
665niter_atsec = 1                        ! Update counter of iterations
666
667zeqn_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
674IF(zeqn_1 > 0._wp) THEN
675   zh_min = zh_1
676ELSEIF(zeqn_1 < 0._wp) THEN
677   zh_max = zh_1
678ELSE ! zh_1 is the root
679   solve_at_general_sec = zh_1
680   IF(PRESENT(p_val)) p_val = zeqn_1
681ENDIF
682
683IF(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
692ELSE
693   zeqn_absmin = ABS(zeqn_1)
694ENDIF
695
696! Pre-calculate the first secant iterate (this is the second iterate)
697niter_atsec = 2
698
699zh_delta = -zeqn_1/((zeqn_2-zeqn_1)/(zh_2 - zh_1))
700zh = 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)
704IF (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)
707ENDIF
708
709IF (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)
712ENDIF
713
714DO
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
801ENDDO
802
803SOLVE_AT_GENERAL_SEC = zh
804
805IF(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
814ENDIF
815
816IF (lhook) CALL dr_hook(RoutineName,zhook_out ,zhook_handle)
817RETURN
818IF (lhook) CALL dr_hook(RoutineName,zhook_out,zhook_handle)
819END FUNCTION SOLVE_AT_GENERAL_SEC
820
821!===============================================================================
822
823FUNCTION 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
831USE mocsy_singledouble
832USE yomhook, ONLY: lhook, dr_hook
833USE parkind1, ONLY: jprb, jpim
834
835IMPLICIT NONE
836REAL(KIND=wp) :: SOLVE_AT_FAST
837
838! Argument variables
839REAL(KIND=wp), INTENT(IN)            :: p_alktot
840REAL(KIND=wp), INTENT(IN)            :: p_dictot
841REAL(KIND=wp), INTENT(IN)            :: p_bortot
842REAL(KIND=wp), INTENT(IN)            :: p_po4tot
843REAL(KIND=wp), INTENT(IN)            :: p_siltot
844!REAL(KIND=wp), INTENT(IN)            :: p_nh4tot
845!REAL(KIND=wp), INTENT(IN)            :: p_h2stot
846REAL(KIND=wp), INTENT(IN)            :: p_so4tot
847REAL(KIND=wp), INTENT(IN)            :: p_flutot
848REAL(KIND=wp), INTENT(IN)            :: K0, K1, K2, Kb, Kw, Ks, Kf
849REAL(KIND=wp), INTENT(IN)            :: K1p, K2p, K3p, Ksi
850REAL(KIND=wp), INTENT(IN), OPTIONAL  :: p_hini
851REAL(KIND=wp), INTENT(OUT), OPTIONAL :: p_val
852
853! Local variables
854REAL(KIND=wp)  ::  zh_ini, zh, zh_prev, zh_lnfactor
855REAL(KIND=wp)  ::  zhdelta
856REAL(KIND=wp)  ::  zeqn, zdeqndh
857!REAL(KIND=wp)  ::  aphscale
858LOGICAL        :: l_exitnow
859REAL(KIND=wp), PARAMETER :: pz_exp_threshold = 1.0_wp
860INTEGER(KIND=jpim), PARAMETER :: zhook_in = 0
861INTEGER(KIND=jpim), PARAMETER :: zhook_out = 1
862REAL(KIND=jprb)               :: zhook_handle
863
864CHARACTER(LEN=*), PARAMETER :: RoutineName='SOLVE_AT_FAST'
865
866IF (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
872IF(PRESENT(p_hini)) THEN
873   zh_ini = p_hini
874ELSE
875   CALL AHINI_FOR_AT(p_alktot, p_dictot, p_bortot, K1, K2, Kb, zh_ini)
876ENDIF
877zh = zh_ini
878
879niter_atfast    = 0                 ! Reset counters of iterations
880DO
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
921ENDDO
922
923SOLVE_AT_FAST = zh
924
925IF(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
934ENDIF
935
936IF (lhook) CALL dr_hook(RoutineName,zhook_out ,zhook_handle)
937RETURN
938IF (lhook) CALL dr_hook(RoutineName,zhook_out,zhook_handle)
939END FUNCTION solve_at_fast
940!===============================================================================
941
942END MODULE mocsy_phsolvers
Note: See TracBrowser for help on using the repository browser.