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.
sedchem.F90 in NEMO/branches/2021/dev_r14383_PISCES_NEWDEV_PISCO/src/TOP/PISCES/SED – NEMO

source: NEMO/branches/2021/dev_r14383_PISCES_NEWDEV_PISCO/src/TOP/PISCES/SED/sedchem.F90 @ 15075

Last change on this file since 15075 was 15075, checked in by aumont, 3 years ago

major update of the sediment module
  • Property svn:keywords set to Id
File size: 33.9 KB
Line 
1MODULE sedchem
2
3   !!======================================================================
4   !!                        ***  Module sedchem  ***
5   !! sediment :   Variable for chemistry of the CO2 cycle
6   !!======================================================================
7   !!   modules used
8   USE sed     ! sediment global variable
9   USE sedarr
10   USE eosbn2, ONLY : neos
11   USE lib_mpp         ! distribued memory computing library
12
13   IMPLICIT NONE
14   PRIVATE
15
16   !! * Accessibility
17   PUBLIC sed_chem
18   PUBLIC ahini_for_at_sed     !
19   PUBLIC solve_at_general_sed !
20
21   ! Maximum number of iterations for each method
22   INTEGER, PARAMETER :: jp_maxniter_atgen    = 20
23   REAL(wp), PARAMETER :: pp_rdel_ah_target = 1.E-4_wp
24
25   !! * Module variables
26   REAL(wp) :: &
27      calcon = 1.03E-2        ! mean calcite concentration [Ca2+] in sea water [mole/kg solution]
28
29   REAL(wp) ::   rgas   = 83.14472      ! universal gas constants
30
31   ! coeff. for density of sea water (Millero & Poisson 1981)
32   REAL(wp), DIMENSION(5)  :: Adsw                       
33   DATA Adsw/8.24493E-1, -4.0899E-3, 7.6438E-5 , -8.246E-7, 5.3875E-9 /
34
35   REAL(wp), DIMENSION(3)  :: Bdsw 
36   DATA Bdsw / -5.72466E-3, 1.0227E-4, -1.6546E-6 /
37
38   REAL(wp)  :: Cdsw = 4.8314E-4
39
40   REAL(wp), DIMENSION(6)  :: Ddsw                   
41   DATA Ddsw / 999.842594 , 6.793952E-2 , -9.095290E-3, 1.001685E-4, -1.120083E-6, 6.536332E-9/
42
43  REAL(wp) :: devk10  = -25.5
44   REAL(wp) :: devk11  = -15.82
45   REAL(wp) :: devk12  = -29.48
46   REAL(wp) :: devk13  = -20.02
47   REAL(wp) :: devk14  = -18.03
48   REAL(wp) :: devk15  = -9.78
49   REAL(wp) :: devk16  = -48.76
50   REAL(wp) :: devk17  = -14.51
51   REAL(wp) :: devk18  = -23.12
52   REAL(wp) :: devk19  = -26.57
53   REAL(wp) :: devk110  = -29.48
54   REAL(wp) :: devk120 = -14.8
55   REAL(wp) :: devk130 = -26.43
56   !
57   REAL(wp) :: devk20  = 0.1271
58   REAL(wp) :: devk21  = -0.0219
59   REAL(wp) :: devk22  = 0.1622
60   REAL(wp) :: devk23  = 0.1119
61   REAL(wp) :: devk24  = 0.0466
62   REAL(wp) :: devk25  = -0.0090
63   REAL(wp) :: devk26  = 0.5304
64   REAL(wp) :: devk27  = 0.1211
65   REAL(wp) :: devk28  = 0.1758
66   REAL(wp) :: devk29  = 0.2020
67   REAL(wp) :: devk210  = 0.1622
68   REAL(wp) :: devk220 = 0.0020
69   REAL(wp) :: devk230 = 0.0889
70   !
71   REAL(wp) :: devk30  = 0.
72   REAL(wp) :: devk31  = 0.
73   REAL(wp) :: devk32  = 2.608E-3
74   REAL(wp) :: devk33  = -1.409e-3
75   REAL(wp) :: devk34  = 0.316e-3
76   REAL(wp) :: devk35  = -0.942e-3
77   REAL(wp) :: devk36  = 0.
78   REAL(wp) :: devk37  = -0.321e-3
79   REAL(wp) :: devk38  = -2.647e-3
80   REAL(wp) :: devk39  = -3.042e-3
81   REAL(wp) :: devk310  = -2.6080e-3
82   REAL(wp) :: devk320 = -0.4E-3
83   REAL(wp) :: devk330 = -0.905E-3
84   !
85   REAL(wp) :: devk40  = -3.08E-3
86   REAL(wp) :: devk41  = 1.13E-3
87   REAL(wp) :: devk42  = -2.84E-3
88   REAL(wp) :: devk43  = -5.13E-3
89   REAL(wp) :: devk44  = -4.53e-3
90   REAL(wp) :: devk45  = -3.91e-3
91   REAL(wp) :: devk46  = -11.76e-3
92   REAL(wp) :: devk47  = -2.67e-3
93   REAL(wp) :: devk48  = -5.15e-3
94   REAL(wp) :: devk49  = -4.08e-3
95   REAL(wp) :: devk410  = -2.84e-3
96   REAL(wp) :: devk420 = 2.89E-3
97   REAL(wp) :: devk430 = -5.03E-3
98   !
99   REAL(wp) :: devk50  = 0.0877E-3
100   REAL(wp) :: devk51  = -0.1475E-3
101   REAL(wp) :: devk52  = 0.
102   REAL(wp) :: devk53  = 0.0794E-3
103   REAL(wp) :: devk54  = 0.09e-3
104   REAL(wp) :: devk55  = 0.054e-3
105   REAL(wp) :: devk56  = 0.3692E-3
106   REAL(wp) :: devk57  = 0.0427e-3
107   REAL(wp) :: devk58  = 0.09e-3
108   REAL(wp) :: devk59  = 0.0714e-3
109   REAL(wp) :: devk510  = 0.0
110   REAL(wp) :: devk520 = 0.054E-3
111   REAL(wp) :: devk530 = 0.0814E-3
112
113   !! * Substitutions
114#  include "do_loop_substitute.h90"
115
116   !! $Id$
117CONTAINS
118
119   SUBROUTINE sed_chem( kt )
120      !!----------------------------------------------------------------------
121      !!                   ***  ROUTINE sed_chem  ***
122      !!
123      !! ** Purpose :   set chemical constants
124      !!
125      !!   History :
126      !!        !  04-10  (N. Emprin, M. Gehlen )  Original code
127      !!        !  06-04  (C. Ethe)  Re-organization
128      !!----------------------------------------------------------------------
129      !!* Arguments
130      INTEGER, INTENT(in) :: kt                     ! time step
131
132      INTEGER  :: ji, jj, ikt
133      REAL(wp) :: ztc, ztc2
134      REAL(wp) :: zsal, zsal15 
135      REAL(wp) :: zdens0, zaw, zbw, zcw   
136      REAL(wp), DIMENSION(jpi,jpj,15) ::   zchem_data
137      !!----------------------------------------------------------------------
138
139
140      IF( ln_timing )  CALL timing_start('sed_chem')
141
142      IF (lwp) WRITE(numsed,*) ' Getting Chemical constants from tracer model at time kt = ', kt
143      IF (lwp) WRITE(numsed,*) ' '
144
145      ! reading variables
146      zchem_data(:,:,:) = rtrn
147
148      IF (ln_sediment_offline) THEN
149         CALL sed_chem_cst
150      ELSE
151         DO_2D( 1, 1, 1, 1 )
152            ikt = mbkt(ji,jj) 
153            IF ( tmask(ji,jj,ikt) == 1 ) THEN
154               zchem_data(ji,jj,1) = ak13  (ji,jj,ikt)
155               zchem_data(ji,jj,2) = ak23  (ji,jj,ikt)
156               zchem_data(ji,jj,3) = akb3  (ji,jj,ikt)
157               zchem_data(ji,jj,4) = akw3  (ji,jj,ikt)
158               zchem_data(ji,jj,5) = aksp  (ji,jj,ikt)
159               zchem_data(ji,jj,6) = borat (ji,jj,ikt)
160               zchem_data(ji,jj,7) = ak1p3 (ji,jj,ikt)
161               zchem_data(ji,jj,8) = ak2p3 (ji,jj,ikt)
162               zchem_data(ji,jj,9) = ak3p3 (ji,jj,ikt)
163               zchem_data(ji,jj,10)= aksi3 (ji,jj,ikt)
164               zchem_data(ji,jj,11)= sio3eq(ji,jj,ikt)
165               zchem_data(ji,jj,12)= aks3  (ji,jj,ikt)
166               zchem_data(ji,jj,13)= akf3  (ji,jj,ikt)
167               zchem_data(ji,jj,14)= sulfat(ji,jj,ikt)
168               zchem_data(ji,jj,15)= fluorid(ji,jj,ikt)
169            ENDIF
170         END_2D
171
172         CALL pack_arr ( jpoce, ak1s  (1:jpoce), zchem_data(1:jpi,1:jpj,1) , iarroce(1:jpoce) )
173         CALL pack_arr ( jpoce, ak2s  (1:jpoce), zchem_data(1:jpi,1:jpj,2) , iarroce(1:jpoce) )
174         CALL pack_arr ( jpoce, akbs  (1:jpoce), zchem_data(1:jpi,1:jpj,3) , iarroce(1:jpoce) )
175         CALL pack_arr ( jpoce, akws  (1:jpoce), zchem_data(1:jpi,1:jpj,4) , iarroce(1:jpoce) )
176         CALL pack_arr ( jpoce, aksps (1:jpoce), zchem_data(1:jpi,1:jpj,5) , iarroce(1:jpoce) )
177         CALL pack_arr ( jpoce, borats(1:jpoce), zchem_data(1:jpi,1:jpj,6) , iarroce(1:jpoce) )
178         CALL pack_arr ( jpoce, ak1ps (1:jpoce), zchem_data(1:jpi,1:jpj,7) , iarroce(1:jpoce) )
179         CALL pack_arr ( jpoce, ak2ps (1:jpoce), zchem_data(1:jpi,1:jpj,8) , iarroce(1:jpoce) )
180         CALL pack_arr ( jpoce, ak3ps (1:jpoce), zchem_data(1:jpi,1:jpj,9) , iarroce(1:jpoce) )
181         CALL pack_arr ( jpoce, aksis (1:jpoce), zchem_data(1:jpi,1:jpj,10), iarroce(1:jpoce) )
182         CALL pack_arr ( jpoce, sieqs (1:jpoce), zchem_data(1:jpi,1:jpj,11), iarroce(1:jpoce) )
183         CALL pack_arr ( jpoce, aks3s (1:jpoce), zchem_data(1:jpi,1:jpj,12), iarroce(1:jpoce) )
184         CALL pack_arr ( jpoce, akf3s (1:jpoce), zchem_data(1:jpi,1:jpj,13), iarroce(1:jpoce) )
185         CALL pack_arr ( jpoce, sulfats(1:jpoce), zchem_data(1:jpi,1:jpj,14), iarroce(1:jpoce) )
186         CALL pack_arr ( jpoce, fluorids(1:jpoce), zchem_data(1:jpi,1:jpj,15), iarroce(1:jpoce) )
187      ENDIF
188
189      DO ji = 1, jpoce
190         ztc     = temp(ji)
191         ztc2    = ztc * ztc
192         zsal    = salt(ji)
193         zsal15  = SQRT( zsal ) * zsal
194
195         ! Density of Sea Water - F(temp,sal) [kg/m3]
196         zdens0 =  Ddsw(1) + Ddsw(2) * ztc + Ddsw(3) * ztc2 &
197                  + Ddsw(4) * ztc * ztc2 + Ddsw(5) * ztc2 * ztc2 &
198                  + Ddsw(6) * ztc * ztc2 * ztc2
199         zaw =  Adsw(1) + Adsw(2) * ztc + Adsw(3)* ztc2 + Adsw(4) * ztc * ztc2 &
200              + Adsw(5) * ztc2 * ztc2
201         zbw =  Bdsw(1) + Bdsw(2) * ztc + Bdsw(3) * ztc2
202         zcw =  Cdsw
203         densSW(ji) = zdens0 + zaw * zsal + zbw * zsal15 + zcw * zsal * zsal
204         densSW(ji) = densSW(ji) * 1E-3   ! to get dens in [kg/l]
205
206         ak12s  (ji) = ak1s (ji) * ak2s (ji)
207         ak12ps (ji) = ak1ps(ji) * ak2ps(ji)
208         ak123ps(ji) = ak1ps(ji) * ak2ps(ji) * ak3ps(ji)
209
210         calcon2(ji) = 0.01028 * ( salt(ji) / 35. ) * densSW(ji)
211      ENDDO
212       
213      IF( ln_timing )  CALL timing_stop('sed_chem')
214
215   END SUBROUTINE sed_chem
216
217   SUBROUTINE ahini_for_at_sed(p_hini)
218      !!---------------------------------------------------------------------
219      !!                     ***  ROUTINE ahini_for_at  ***
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
223      !! polynomial) around the local minimum, if it exists.
224      !! Returns * 1E-03_wp if p_alkcb <= 0
225      !!         * 1E-10_wp if p_alkcb >= 2*p_dictot + p_bortot
226      !!         * 1E-07_wp if 0 < p_alkcb < 2*p_dictot + p_bortot
227      !!                    and the 2nd order approximation does not have
228      !!                    a solution
229      !!---------------------------------------------------------------------
230      REAL(wp), DIMENSION(jpoce,jpksed), INTENT(OUT)  ::  p_hini
231      INTEGER  ::   ji, jk
232      REAL(wp)  ::  zca1, zba1
233      REAL(wp)  ::  zd, zsqrtd, zhmin
234      REAL(wp)  ::  za2, za1, za0
235      REAL(wp)  ::  p_dictot, p_bortot, p_alkcb
236
237      IF( ln_timing )  CALL timing_start('ahini_for_at_sed')
238      !
239      DO jk = 1, jpksed
240         DO ji = 1, jpoce
241            p_alkcb  = pwcp(ji,jk,jwalk) / densSW(ji)
242            p_dictot = pwcp(ji,jk,jwdic) / densSW(ji) 
243            p_bortot = borats(ji)
244            IF (p_alkcb <= 0.) THEN
245                p_hini(ji,jk) = 1.e-3
246            ELSEIF (p_alkcb >= (2.*p_dictot + p_bortot)) THEN
247                p_hini(ji,jk) = 1.e-10_wp
248            ELSE
249                zca1 = p_dictot/( p_alkcb + rtrn )
250                zba1 = p_bortot/ (p_alkcb + rtrn )
251           ! Coefficients of the cubic polynomial
252                za2 = aKbs(ji)*(1. - zba1) + ak1s(ji)*(1.-zca1)
253                za1 = ak1s(ji)*akbs(ji)*(1. - zba1 - zca1)    &
254                &     + ak1s(ji)*ak2s(ji)*(1. - (zca1+zca1))
255                za0 = ak1s(ji)*ak2s(ji)*akbs(ji)*(1. - zba1 - (zca1+zca1))
256                                        ! Taylor expansion around the minimum
257                zd = za2*za2 - 3.*za1   ! Discriminant of the quadratic equation
258                                        ! for the minimum close to the root
259
260                IF(zd > 0.) THEN        ! If the discriminant is positive
261                  zsqrtd = SQRT(zd)
262                  IF(za2 < 0) THEN
263                    zhmin = (-za2 + zsqrtd)/3.
264                  ELSE
265                    zhmin = -za1/(za2 + zsqrtd)
266                  ENDIF
267                  p_hini(ji,jk) = zhmin + SQRT(-(za0 + zhmin*(za1 + zhmin*(za2 + zhmin)))/zsqrtd)
268                ELSE
269                  p_hini(ji,jk) = 1.e-7
270                ENDIF
271             !
272            ENDIF
273         END DO
274      END DO
275      !
276      IF( ln_timing )  CALL timing_stop('ahini_for_at_sed')
277      !
278   END SUBROUTINE ahini_for_at_sed
279
280   !===============================================================================
281   SUBROUTINE anw_infsup_sed( p_alknw_inf, p_alknw_sup )
282
283   ! Subroutine returns the lower and upper bounds of "non-water-selfionization"
284   ! contributions to total alkalinity (the infimum and the supremum), i.e
285   ! inf(TA - [OH-] + [H+]) and sup(TA - [OH-] + [H+])
286
287   ! Argument variables
288   INTEGER :: jk
289   REAL(wp), DIMENSION(jpoce,jpksed), INTENT(OUT) :: p_alknw_inf
290   REAL(wp), DIMENSION(jpoce,jpksed), INTENT(OUT) :: p_alknw_sup
291
292   DO jk = 1, jpksed
293      p_alknw_inf(:,jk) =  -pwcp(:,jk,jwpo4) / densSW(:)
294      p_alknw_sup(:,jk) =   (2. * pwcp(:,jk,jwdic) + 2. * pwcp(:,jk,jwpo4) + pwcp(:,jk,jwsil) )    &
295   &                        / densSW(:) + borats(:) 
296   END DO
297
298   END SUBROUTINE anw_infsup_sed
299
300
301   SUBROUTINE solve_at_general_sed( p_hini, zhi )
302
303   ! Universal pH solver that converges from any given initial value,
304   ! determines upper an lower bounds for the solution if required
305
306   ! Argument variables
307   !--------------------
308   REAL(wp), DIMENSION(jpoce,jpksed), INTENT(IN)   :: p_hini
309   REAL(wp), DIMENSION(jpoce,jpksed), INTENT(OUT)  :: zhi
310
311   ! Local variables
312   !-----------------
313   INTEGER   ::  ji, jk, jn
314   REAL(wp)  ::  zh_ini, zh, zh_prev, zh_lnfactor
315   REAL(wp)  ::  zdelta, zh_delta
316   REAL(wp)  ::  zeqn, zdeqndh, zalka
317   REAL(wp)  ::  aphscale
318   REAL(wp)  ::  znumer_dic, zdnumer_dic, zdenom_dic, zalk_dic, zdalk_dic
319   REAL(wp)  ::  znumer_bor, zdnumer_bor, zdenom_bor, zalk_bor, zdalk_bor
320   REAL(wp)  ::  znumer_po4, zdnumer_po4, zdenom_po4, zalk_po4, zdalk_po4
321   REAL(wp)  ::  znumer_sil, zdnumer_sil, zdenom_sil, zalk_sil, zdalk_sil
322   REAL(wp)  ::  znumer_so4, zdnumer_so4, zdenom_so4, zalk_so4, zdalk_so4
323   REAL(wp)  ::  znumer_flu, zdnumer_flu, zdenom_flu, zalk_flu, zdalk_flu
324   REAL(wp)  ::  znumer_h2s, zdnumer_h2s, zdenom_h2s, zalk_h2s, zdalk_h2s
325   REAL(wp)  ::  znumer_nh3, zdnumer_nh3, zdenom_nh3, zalk_nh3, zdalk_nh3
326   REAL(wp)  ::  zalk_wat, zdalk_wat
327   REAL(wp)  ::  zfact, p_alktot, zdic, zbot, zpt, zst, zft, zsit, zh2s, znh3
328   LOGICAL   ::  l_exitnow
329   REAL(wp), PARAMETER :: pz_exp_threshold = 1.0
330   REAL(wp), DIMENSION(jpoce,jpksed) :: zalknw_inf, zalknw_sup, rmask, zh_min, zh_max, zeqn_absmin
331
332   IF( ln_timing )  CALL timing_start('solve_at_general_sed')
333      !  Allocate temporary workspace
334   CALL anw_infsup_sed( zalknw_inf, zalknw_sup )
335
336   rmask(:,:) = 1.0
337   zhi(:,:)   = 0.
338
339   ! TOTAL H+ scale: conversion factor for Htot = aphscale * Hfree
340   DO jk = 1, jpksed
341      DO ji = 1, jpoce
342         IF (rmask(ji,jk) == 1.) THEN
343            p_alktot = pwcp(ji,jk,jwalk) / densSW(ji)
344            aphscale = 1. + sulfats(ji)/aks3s(ji)
345            zh_ini = p_hini(ji,jk)
346
347            zdelta = (p_alktot-zalknw_inf(ji,jk))**2 + 4.*akws(ji) / aphscale
348
349            IF(p_alktot >= zalknw_inf(ji,jk)) THEN
350               zh_min(ji,jk) = 2.*akws(ji) /( p_alktot-zalknw_inf(ji,jk) + SQRT(zdelta) )
351            ELSE
352               zh_min(ji,jk) = aphscale * (-(p_alktot-zalknw_inf(ji,jk)) + SQRT(zdelta) ) / 2.
353            ENDIF
354
355            zdelta = (p_alktot-zalknw_sup(ji,jk))**2 + 4.*akws(ji) / aphscale
356
357            IF(p_alktot <= zalknw_sup(ji,jk)) THEN
358               zh_max(ji,jk) = aphscale * (-(p_alktot-zalknw_sup(ji,jk)) + SQRT(zdelta) ) / 2.
359            ELSE
360               zh_max(ji,jk) = 2.*akws(ji) /( p_alktot-zalknw_sup(ji,jk) + SQRT(zdelta) )
361            ENDIF
362
363            zhi(ji,jk) = MAX(MIN(zh_max(ji,jk), zh_ini), zh_min(ji,jk))
364         ENDIF
365      END DO
366   END DO
367
368   zeqn_absmin(:,:) = HUGE(1._wp)
369
370   DO jn = 1, jp_maxniter_atgen
371   DO jk = 1, jpksed
372      DO ji = 1, jpoce
373         IF (rmask(ji,jk) == 1.) THEN
374
375            p_alktot = pwcp(ji,jk,jwalk) / densSW(ji)
376            zdic = pwcp(ji,jk,jwdic) / densSW(ji)
377            zbot = borats(ji)
378            zpt  = pwcp(ji,jk,jwpo4) / densSW(ji)
379            zsit = pwcp(ji,jk,jwsil) / densSW(ji)
380            zst  = pwcp(ji,jk,jwso4) / densSW(ji)
381            zft  = fluorids(ji)
382            zh2s = pwcp(ji,jk,jwh2s) / densSW(ji)
383            znh3 = pwcp(ji,jk,jwnh4) / densSW(ji)
384            aphscale = 1. + sulfats(ji)/aks3s(ji)
385            zh = zhi(ji,jk)
386            zh_prev = zh
387
388            ! H2CO3 - HCO3 - CO3 : n=2, m=0
389            znumer_dic = 2.*ak1s(ji)*ak2s(ji) + zh*ak1s(ji)
390            zdenom_dic = ak1s(ji)*ak2s(ji) + zh*(ak1s(ji) + zh)
391            zalk_dic   = zdic * (znumer_dic/zdenom_dic)
392            zdnumer_dic = ak1s(ji)*ak1s(ji)*ak2s(ji) + zh     &
393                          *(4.*ak1s(ji)*ak2s(ji) + zh*ak1s(ji))
394            zdalk_dic   = -zdic*(zdnumer_dic/zdenom_dic**2)
395
396
397            ! B(OH)3 - B(OH)4 : n=1, m=0
398            znumer_bor  = akbs(ji)
399            zdenom_bor  = akbs(ji) + zh
400            zalk_bor    = zbot * (znumer_bor/zdenom_bor)
401            zdnumer_bor = akbs(ji)
402            zdalk_bor   = -zbot*(zdnumer_bor/zdenom_bor**2)
403
404
405            ! H3PO4 - H2PO4 - HPO4 - PO4 : n=3, m=1
406            znumer_po4 = 3.*ak1ps(ji)*ak2ps(ji)*ak3ps(ji)  &
407            &            + zh*(2.*ak1ps(ji)*ak2ps(ji) + zh* ak1ps(ji))
408            zdenom_po4 = ak1ps(ji)*ak2ps(ji)*ak3ps(ji)     &
409            &            + zh*( ak1ps(ji)*ak2ps(ji) + zh*(ak1ps(ji) + zh))
410            zalk_po4   = zpt * (znumer_po4/zdenom_po4 - 1.) ! Zero level of H3PO4 = 1
411            zdnumer_po4 = ak1ps(ji)*ak2ps(ji)*ak1ps(ji)*ak2ps(ji)*ak3ps(ji)  &
412            &             + zh*(4.*ak1ps(ji)*ak1ps(ji)*ak2ps(ji)*ak3ps(ji)         &
413            &             + zh*(9.*ak1ps(ji)*ak2ps(ji)*ak3ps(ji)                         &
414            &             + ak1ps(ji)*ak1ps(ji)*ak2ps(ji)                                &
415            &             + zh*(4.*ak1ps(ji)*ak2ps(ji) + zh * ak1ps(ji) ) ) )
416            zdalk_po4   = -zpt * (zdnumer_po4/zdenom_po4**2)
417
418            ! H4SiO4 - H3SiO4 : n=1, m=0
419            znumer_sil = aksis(ji)
420            zdenom_sil = aksis(ji) + zh
421            zalk_sil   = zsit * (znumer_sil/zdenom_sil)
422            zdnumer_sil = aksis(ji)
423            zdalk_sil   = -zsit * (zdnumer_sil/zdenom_sil**2)
424
425            ! H2S - HS : n=1i, m=1
426            znumer_h2s  = akh2s(ji)
427            zdenom_h2s  = akh2s(ji) + zh
428            zalk_h2s    = zh2s * (znumer_h2s/zdenom_h2s)
429            zdnumer_h2s = akh2s(ji)
430            zdalk_h2s   = -zh2s * (zdnumer_h2s/zdenom_h2s**2)
431
432            ! NH4 - NH3 : n=1i, m=1
433            znumer_nh3  = aknh3(ji)
434            zdenom_nh3  = aknh3(ji) + zh
435            zalk_nh3    = znh3 * (znumer_nh3/zdenom_nh3)
436            zdnumer_nh3 = aknh3(ji)
437            zdalk_nh3   = -znh3 * (zdnumer_nh3/zdenom_nh3**2)
438
439            ! HSO4 - SO4 : n=1, m=1
440            aphscale = 1.0 + zst/aks3s(ji)
441            znumer_so4 = aks3s(ji) * aphscale
442            zdenom_so4 = aks3s(ji) * aphscale + zh
443            zalk_so4   = zst * (znumer_so4/zdenom_so4 - 1.)
444            zdnumer_so4 = aks3s(ji) * aphscale
445            zdalk_so4   = -zst * (zdnumer_so4/zdenom_so4**2)
446
447            ! HF - F : n=1, m=1
448            znumer_flu =  akf3s(ji)
449            zdenom_flu =  akf3s(ji) + zh
450            zalk_flu   =  zft * (znumer_flu/zdenom_flu - 1.)
451            zdnumer_flu = akf3s(ji)
452            zdalk_flu   = -zft * (zdnumer_flu/zdenom_flu**2)
453
454            ! H2O - OH
455            zalk_wat   = akws(ji)/zh - zh/aphscale
456            zdalk_wat  = -akws(ji)/zh**2 - 1./aphscale
457
458            ! CALCULATE [ALK]([CO3--], [HCO3-])
459            zeqn = zalk_dic + zalk_bor + zalk_po4 + zalk_sil   &
460            &      + zalk_so4 + zalk_flu                       &
461            &      + zalk_wat + zalk_h2s + zalk_nh3 - p_alktot
462
463            zalka = p_alktot - (zalk_bor + zalk_po4 + zalk_sil   &
464            &       + zalk_so4 + zalk_flu + zalk_wat + zalk_h2s + zalk_nh3)
465
466            zdeqndh = zdalk_dic + zdalk_bor + zdalk_po4 + zdalk_sil &
467            &         + zdalk_so4 + zdalk_flu + zdalk_wat + zdalk_h2s + zdalk_nh3
468
469            ! Adapt bracketing interval
470            IF(zeqn > 0._wp) THEN
471               zh_min(ji,jk) = zh_prev
472            ELSEIF(zeqn < 0._wp) THEN
473               zh_max(ji,jk) = zh_prev
474            ENDIF
475
476            IF(ABS(zeqn) >= 0.5_wp*zeqn_absmin(ji,jk)) THEN
477            ! if the function evaluation at the current point is
478            ! not decreasing faster than with a bisection step (at least linearly)
479            ! in absolute value take one bisection step on [ph_min, ph_max]
480            ! ph_new = (ph_min + ph_max)/2d0
481            !
482            ! In terms of [H]_new:
483            ! [H]_new = 10**(-ph_new)
484            !         = 10**(-(ph_min + ph_max)/2d0)
485            !         = SQRT(10**(-(ph_min + phmax)))
486            !         = SQRT(zh_max * zh_min)
487               zh = SQRT(zh_max(ji,jk) * zh_min(ji,jk))
488               zh_lnfactor = (zh - zh_prev)/zh_prev ! Required to test convergence below
489            ELSE
490            ! dzeqn/dpH = dzeqn/d[H] * d[H]/dpH
491            !           = -zdeqndh * LOG(10) * [H]
492            ! \Delta pH = -zeqn/(zdeqndh*d[H]/dpH) = zeqn/(zdeqndh*[H]*LOG(10))
493            !
494            ! pH_new = pH_old + \deltapH
495            !
496            ! [H]_new = 10**(-pH_new)
497            !         = 10**(-pH_old - \Delta pH)
498            !         = [H]_old * 10**(-zeqn/(zdeqndh*[H]_old*LOG(10)))
499            !         = [H]_old * EXP(-LOG(10)*zeqn/(zdeqndh*[H]_old*LOG(10)))
500            !         = [H]_old * EXP(-zeqn/(zdeqndh*[H]_old))
501
502               zh_lnfactor = -zeqn/(zdeqndh*zh_prev)
503
504               IF(ABS(zh_lnfactor) > pz_exp_threshold) THEN
505                  zh          = zh_prev*EXP(zh_lnfactor)
506               ELSE
507                  zh_delta    = zh_lnfactor*zh_prev
508                  zh          = zh_prev + zh_delta
509               ENDIF
510
511               IF( zh < zh_min(ji,jk) ) THEN
512               ! if [H]_new < [H]_min
513               ! i.e., if ph_new > ph_max then
514               ! take one bisection step on [ph_prev, ph_max]
515               ! ph_new = (ph_prev + ph_max)/2d0
516               ! In terms of [H]_new:
517               ! [H]_new = 10**(-ph_new)
518               !         = 10**(-(ph_prev + ph_max)/2d0)
519               !         = SQRT(10**(-(ph_prev + phmax)))
520               !         = SQRT([H]_old*10**(-ph_max))
521               !         = SQRT([H]_old * zh_min)
522                  zh                = SQRT(zh_prev * zh_min(ji,jk))
523                  zh_lnfactor       = (zh - zh_prev)/zh_prev ! Required to test convergence below
524               ENDIF
525
526               IF( zh > zh_max(ji,jk) ) THEN
527               ! if [H]_new > [H]_max
528               ! i.e., if ph_new < ph_min, then
529               ! take one bisection step on [ph_min, ph_prev]
530               ! ph_new = (ph_prev + ph_min)/2d0
531               ! In terms of [H]_new:
532               ! [H]_new = 10**(-ph_new)
533               !         = 10**(-(ph_prev + ph_min)/2d0)
534               !         = SQRT(10**(-(ph_prev + ph_min)))
535               !         = SQRT([H]_old*10**(-ph_min))
536               !         = SQRT([H]_old * zhmax)
537                  zh                = SQRT(zh_prev * zh_max(ji,jk))
538                  zh_lnfactor       = (zh - zh_prev)/zh_prev ! Required to test convergence below
539               ENDIF
540            ENDIF
541
542            zeqn_absmin(ji,jk) = MIN( ABS(zeqn), zeqn_absmin(ji,jk))
543
544            ! Stop iterations once |\delta{[H]}/[H]| < rdel
545            ! <=> |(zh - zh_prev)/zh_prev| = |EXP(-zeqn/(zdeqndh*zh_prev)) -1| < rdel
546            ! |EXP(-zeqn/(zdeqndh*zh_prev)) -1| ~ |zeqn/(zdeqndh*zh_prev)|
547            ! Alternatively:
548            ! |\Delta pH| = |zeqn/(zdeqndh*zh_prev*LOG(10))|
549            !             ~ 1/LOG(10) * |\Delta [H]|/[H]
550            !             < 1/LOG(10) * rdel
551
552            ! Hence |zeqn/(zdeqndh*zh)| < rdel
553
554            ! rdel <-- pp_rdel_ah_target
555            l_exitnow = (ABS(zh_lnfactor) < pp_rdel_ah_target)
556
557            IF(l_exitnow) THEN
558               rmask(ji,jk) = 0.
559            ENDIF
560
561            zhi(ji,jk) =  zh
562
563            IF(jn >= jp_maxniter_atgen) THEN
564               zhi(ji,jk) = -1._wp
565            ENDIF
566
567         ENDIF
568      END DO
569   END DO
570   END DO
571   !
572   IF( ln_timing )  CALL timing_stop('solve_at_general_sed')
573
574   END SUBROUTINE solve_at_general_sed
575
576   SUBROUTINE sed_chem_cst
577      !!---------------------------------------------------------------------
578      !!                     ***  ROUTINE sed_chem_cst  ***
579      !!
580      !! ** Purpose :   Sea water chemistry computed following MOCSY protocol
581      !!                Computation is done at the bottom of the ocean only
582      !!
583      !! ** Method  : - ...
584      !!---------------------------------------------------------------------
585      INTEGER  ::   ji
586      REAL(wp), DIMENSION(jpoce) :: saltprac, temps
587      REAL(wp) ::   ztkel, ztkel1, zt , zsal  , zsal2 , zbuf1 , zbuf2
588      REAL(wp) ::   ztgg , ztgg2, ztgg3 , ztgg4 , ztgg5
589      REAL(wp) ::   zpres, ztc  , zcl   , zcpexp, zoxy  , zcpexp2
590      REAL(wp) ::   zsqrt, ztr  , zlogt , zcek1, zc1, zplat
591      REAL(wp) ::   zis  , zis2 , zsal15, zisqrt, za1, za2
592      REAL(wp) ::   zckb , zck1 , zck2  , zckw  , zak1 , zak2  , zakb , zaksp0, zakw
593      REAL(wp) ::   zck1p, zck2p, zck3p, zcksi, zak1p, zak2p, zak3p, zaksi
594      REAL(wp) ::   zst  , zft  , zcks  , zckhs, zckf  , zaksp1, zcknh3
595      REAL(wp) ::   total2free, free2SWS, total2SWS, SWS2total
596      !!---------------------------------------------------------------------
597      !
598      IF( ln_timing )   CALL timing_start('sed_chem_cst')
599      !
600      ! Computation of chemical constants require practical salinity
601      ! Thus, when TEOS08 is used, absolute salinity is converted to
602      ! practical salinity
603      ! -------------------------------------------------------------
604      IF (neos == -1) THEN
605         saltprac(:) = salt(:) * 35.0 / 35.16504
606      ELSE
607         saltprac(:) = salt(:)
608      ENDIF
609
610      !
611      ! Computations of chemical constants require in situ temperature
612      ! Here a quite simple formulation is used to convert
613      ! potential temperature to in situ temperature. The errors is less than
614      ! 0.04°C relative to an exact computation
615      ! ---------------------------------------------------------------------
616         DO ji = 1, jpoce
617            zpres = zkbot(ji) / 1000.
618            za1 = 0.04 * ( 1.0 + 0.185 * temp(ji) + 0.035 * (saltprac(ji) - 35.0) )
619            za2 = 0.0075 * ( 1.0 - temp(ji) / 30.0 )
620            temps(ji) = temp(ji) - za1 * zpres + za2 * zpres**2
621         END DO
622
623      ! CHEMICAL CONSTANTS - DEEP OCEAN
624      ! -------------------------------
625      DO ji = 1, jpoce
626         ! SET PRESSION ACCORDING TO SAUNDER (1980)
627         zc1 = 5.92E-3 
628         zpres = ((1-zc1)-SQRT(((1-zc1)**2)-(8.84E-6*zkbot(ji)))) / 4.42E-6
629         zpres = zpres / 10.0
630
631         ! SET ABSOLUTE TEMPERATURE
632         ztkel   = temps(ji) + 273.15
633         zsal    = saltprac(ji)
634         zsqrt  = SQRT( zsal )
635         zsal15  = zsqrt * zsal
636         zlogt  = LOG( ztkel )
637         ztr    = 1. / ztkel
638         zis    = 19.924 * zsal / ( 1000.- 1.005 * zsal )
639         zis2   = zis * zis
640         zisqrt = SQRT( zis )
641         ztc    = temps(ji)
642
643         ! CHLORINITY (WOOSTER ET AL., 1969)
644         zcl     = zsal / 1.80655
645
646         ! TOTAL SULFATE CONCENTR. [MOLES/kg soln]
647         zst     = 0.14 * zcl /96.062
648
649         ! TOTAL FLUORIDE CONCENTR. [MOLES/kg soln]
650         zft     = 0.000067 * zcl /18.9984
651
652         ! DISSOCIATION CONSTANT FOR SULFATES on free H scale (Dickson 1990)
653         zcks    = EXP(-4276.1 * ztr + 141.328 - 23.093 * zlogt         &
654         &         + (-13856. * ztr + 324.57 - 47.986 * zlogt) * zisqrt &
655         &         + (35474. * ztr - 771.54 + 114.723 * zlogt) * zis    &
656         &         - 2698. * ztr * zis**1.5 + 1776.* ztr * zis2         &
657         &         + LOG(1.0 - 0.001005 * zsal))
658
659         ! DISSOCIATION CONSTANT FOR H2S on free H scale
660         zckhs   = EXP(-13275.3 * ztr + 225.838 - 34.6435 * zlogt       &
661         &         + 0.3449 * zsqrt - 0.0274 * zsal )
662
663         ! DISSOCIATION CONSTANT FOR NH3 on free H scale
664         zcknh3   = EXP(-6285.33 * ztr - 0.25444  + 0.0001635 * ztkel   &
665         &         + (0.46532 - 123.7184 * ztr) * zsqrt                 &
666         &         + (-0.01992 + 3.17556 * ztr) * zsal )
667
668         ! DISSOCIATION CONSTANT FOR FLUORIDES on free H scale (Dickson and Riley 79)
669         zckf    = EXP( 1590.2*ztr - 12.641 + 1.525*zisqrt   &
670         &         + LOG(1.0d0 - 0.001005d0*zsal)            &
671         &         + LOG(1.0d0 + zst/zcks))
672
673         ! DISSOCIATION CONSTANT FOR CARBONATE AND BORATE
674         zckb=  (-8966.90 - 2890.53*zsqrt - 77.942*zsal        &
675         &      + 1.728*zsal15 - 0.0996*zsal*zsal)*ztr         &
676         &      + (148.0248 + 137.1942*zsqrt + 1.62142*zsal)   &
677         &      + (-24.4344 - 25.085*zsqrt - 0.2474*zsal)      &
678         &      * zlogt + 0.053105*zsqrt*ztkel
679
680         ! DISSOCIATION COEFFICIENT FOR CARBONATE ACCORDING TO
681         ! MEHRBACH (1973) REFIT BY MILLERO (1995), seawater scale
682         zck1    = -1.0*(3633.86*ztr - 61.2172 + 9.6777*zlogt  &
683                   - 0.011555*zsal + 0.0001152*zsal*zsal)
684         zck2    = -1.0*(471.78*ztr + 25.9290 - 3.16967*zlogt      &
685                   - 0.01781*zsal + 0.0001122*zsal*zsal)
686
687         ! PKW (H2O) (MILLERO, 1995) from composite data
688         zckw    = -13847.26 * ztr + 148.9652 - 23.6521 * zlogt + ( 118.67 * ztr    &
689                   - 5.977 + 1.0495 * zlogt ) * zsqrt - 0.01615 * zsal
690
691         ! CONSTANTS FOR PHOSPHATE (MILLERO, 1995)
692         zck1p    = -4576.752*ztr + 115.540 - 18.453*zlogt   &
693         &          + (-106.736*ztr + 0.69171) * zsqrt       &
694         &          + (-0.65643*ztr - 0.01844) * zsal
695
696         zck2p    = -8814.715*ztr + 172.1033 - 27.927*zlogt  &
697         &          + (-160.340*ztr + 1.3566)*zsqrt          &
698         &          + (0.37335*ztr - 0.05778)*zsal
699
700         zck3p    = -3070.75*ztr - 18.126                    &
701         &          + (17.27039*ztr + 2.81197) * zsqrt       &
702         &          + (-44.99486*ztr - 0.09984) * zsal
703
704         ! CONSTANT FOR SILICATE, MILLERO (1995)
705         zcksi    = -8904.2*ztr  + 117.400 - 19.334*zlogt   &
706         &          + (-458.79*ztr + 3.5913) * zisqrt       &
707         &          + (188.74*ztr - 1.5998) * zis           &
708         &          + (-12.1652*ztr + 0.07871) * zis2       &
709         &          + LOG(1.0 - 0.001005*zsal)
710
711         ! APPARENT SOLUBILITY PRODUCT K'SP OF CALCITE IN SEAWATER
712         !       (S=27-43, T=2-25 DEG C) at pres =0 (atmos. pressure) (MUCCI 1983)
713         zaksp0  = -171.9065 -0.077993*ztkel + 2839.319*ztr + 71.595*LOG10( ztkel )   &
714         &         + (-0.77712 + 0.00284263*ztkel + 178.34*ztr) * zsqrt  &
715         &         - 0.07711*zsal + 0.0041249*zsal15
716
717         ! CONVERT FROM DIFFERENT PH SCALES
718         total2free  = 1.0/(1.0 + zst/zcks)
719         free2SWS    = 1. + zst/zcks + zft/(zckf*total2free)
720         total2SWS   = total2free * free2SWS
721         SWS2total   = 1.0 / total2SWS
722
723
724         ! K1, K2 OF CARBONIC ACID, KB OF BORIC ACID, KW (H2O) (LIT.?)
725         zak1    = 10**(zck1) * total2SWS
726         zak2    = 10**(zck2) * total2SWS
727         zakb    = EXP( zckb ) * total2SWS
728         zakw    = EXP( zckw )
729         zaksp1  = 10**(zaksp0)
730         zak1p   = exp( zck1p )
731         zak2p   = exp( zck2p )
732         zak3p   = exp( zck3p )
733         zaksi   = exp( zcksi )
734         zckf    = zckf * total2SWS
735
736         ! FORMULA FOR CPEXP AFTER EDMOND & GIESKES (1970)
737         !        (REFERENCE TO CULBERSON & PYTKOQICZ (1968) AS MADE
738         !        IN BROECKER ET AL. (1982) IS INCORRECT; HERE RGAS IS
739         !        TAKEN TENFOLD TO CORRECT FOR THE NOTATION OF pres  IN
740         !        DBAR INSTEAD OF BAR AND THE EXPRESSION FOR CPEXP IS
741         !        MULTIPLIED BY LN(10.) TO ALLOW USE OF EXP-FUNCTION
742         !        WITH BASIS E IN THE FORMULA FOR AKSPP (CF. EDMOND
743         !        & GIESKES (1970), P. 1285-1286 (THE SMALL
744         !        FORMULA ON P. 1286 IS RIGHT AND CONSISTENT WITH THE
745         !        SIGN IN PARTIAL MOLAR VOLUME CHANGE AS SHOWN ON P. 1285))
746         zcpexp  = zpres / (rgas*ztkel)
747         zcpexp2 = zpres * zcpexp
748
749         ! KB OF BORIC ACID, K1,K2 OF CARBONIC ACID PRESSURE
750         !        CORRECTION AFTER CULBERSON AND PYTKOWICZ (1968)
751         !        (CF. BROECKER ET AL., 1982)
752
753         zbuf1  = -     ( devk10 + devk20 * ztc + devk30 * ztc * ztc )
754         zbuf2  = 0.5 * ( devk40 + devk50 * ztc )
755         ak1s(ji) = zak1 * EXP( zbuf1 * zcpexp + zbuf2 * zcpexp2 )
756
757         zbuf1  =     - ( devk11 + devk21 * ztc + devk31 * ztc * ztc )
758         zbuf2  = 0.5 * ( devk41 + devk51 * ztc )
759         ak2s(ji) = zak2 * EXP( zbuf1 * zcpexp + zbuf2 * zcpexp2 )
760
761         zbuf1  =     - ( devk12 + devk22 * ztc + devk32 * ztc * ztc )
762         zbuf2  = 0.5 * ( devk42 + devk52 * ztc )
763         akbs(ji) = zakb * EXP( zbuf1 * zcpexp + zbuf2 * zcpexp2 )
764
765         zbuf1  =     - ( devk13 + devk23 * ztc + devk33 * ztc * ztc )
766         zbuf2  = 0.5 * ( devk43 + devk53 * ztc )
767         akws(ji) = zakw * EXP( zbuf1 * zcpexp + zbuf2 * zcpexp2 )
768
769         zbuf1  =     - ( devk14 + devk24 * ztc + devk34 * ztc * ztc )
770         zbuf2  = 0.5 * ( devk44 + devk54 * ztc )
771         aks3s(ji) = zcks * EXP( zbuf1 * zcpexp + zbuf2 * zcpexp2 )
772
773         zbuf1  =     - ( devk15 + devk25 * ztc + devk35 * ztc * ztc )
774         zbuf2  = 0.5 * ( devk45 + devk55 * ztc )
775         akf3s(ji) = zckf * EXP( zbuf1 * zcpexp + zbuf2 * zcpexp2 )
776
777         zbuf1  =     - ( devk17 + devk27 * ztc + devk37 * ztc * ztc )
778         zbuf2  = 0.5 * ( devk47 + devk57 * ztc )
779         ak1ps(ji) = zak1p * EXP( zbuf1 * zcpexp + zbuf2 * zcpexp2 )
780
781         zbuf1  =     - ( devk18 + devk28 * ztc + devk38 * ztc * ztc )
782         zbuf2  = 0.5 * ( devk48 + devk58 * ztc )
783         ak2ps(ji) = zak2p * EXP( zbuf1 * zcpexp + zbuf2 * zcpexp2 )
784
785         zbuf1  =     - ( devk110 + devk210 * ztc + devk310 * ztc * ztc )
786         zbuf2  = 0.5 * ( devk410 + devk510 * ztc )
787         aksis(ji) = zaksi * EXP( zbuf1 * zcpexp + zbuf2 * zcpexp2 )
788
789         zbuf1  =     - ( devk120 + devk220 * ztc + devk320 * ztc * ztc )
790         zbuf2  = 0.5 * ( devk420 + devk520 * ztc )
791         akh2s(ji) = zckhs * EXP( zbuf1 * zcpexp + zbuf2 * zcpexp2 )
792
793         zbuf1  =     - ( devk130 + devk230 * ztc + devk330 * ztc * ztc )
794         zbuf2  = 0.5 * ( devk430 + devk530 * ztc )
795         aknh3(ji) = zcknh3 * EXP( zbuf1 * zcpexp + zbuf2 * zcpexp2 )
796
797         ! Convert to total scale
798         ak1s(ji)  = ak1s(ji)  * SWS2total
799         ak2s(ji)  = ak2s(ji)  * SWS2total
800         akbs(ji)  = akbs(ji)  * SWS2total
801         akws(ji)  = akws(ji)  * SWS2total
802         ak1ps(ji) = ak1ps(ji) * SWS2total
803         ak2ps(ji) = ak2ps(ji) * SWS2total
804         ak3ps(ji) = ak3ps(ji) * SWS2total
805         aksis(ji) = aksis(ji) * SWS2total
806         akf3s(ji) = akf3s(ji) / total2free
807         akh2s(ji) = akh2s(ji) * SWS2total
808         aknh3(ji) = aknh3(ji) * SWS2total
809
810         ! APPARENT SOLUBILITY PRODUCT K'SP OF CALCITE
811         !        AS FUNCTION OF PRESSURE FOLLOWING MILLERO
812         !        (P. 1285) AND BERNER (1976)
813         zbuf1  =     - ( devk16 + devk26 * ztc + devk36 * ztc * ztc )
814         zbuf2  = 0.5 * ( devk46 + devk56 * ztc )
815         aksps(ji) = zaksp1 * EXP( zbuf1 * zcpexp + zbuf2 * zcpexp2 )
816
817         ! TOTAL F, S, and BORATE CONCENTR. [MOLES/L]
818         borats(ji)   = 0.0002414 * zcl / 10.811
819         sulfats(ji)  = zst
820         fluorids(ji) = zft
821
822         ! Iron and SIO3 saturation concentration from ...
823         sieqs(ji) = EXP(  LOG( 10.) * ( 6.44 - 968. / ztkel )  ) * 1.e-6
824      END DO
825      !
826      IF( ln_timing )  CALL timing_stop('sed_chem_cst')
827      !
828   END SUBROUTINE sed_chem_cst
829
830
831END MODULE sedchem
Note: See TracBrowser for help on using the repository browser.