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/2019/dev_r11943_MERGE_2019/src/TOP/PISCES/SED – NEMO

source: NEMO/branches/2019/dev_r11943_MERGE_2019/src/TOP/PISCES/SED/sedchem.F90 @ 12340

Last change on this file since 12340 was 12340, checked in by acc, 4 years ago

Branch 2019/dev_r11943_MERGE_2019. This commit introduces basic do loop macro
substitution to the 2019 option 1, merge branch. These changes have been SETTE
tested. The only addition is the do_loop_substitute.h90 file in the OCE directory but
the macros defined therein are used throughout the code to replace identifiable, 2D-
and 3D- nested loop opening and closing statements with single-line alternatives. Code
indents are also adjusted accordingly.

The following explanation is taken from comments in the new header file:

This header file contains preprocessor definitions and macros used in the do-loop
substitutions introduced between version 4.0 and 4.2. The primary aim of these macros
is to assist in future applications of tiling to improve performance. This is expected
to be achieved by alternative versions of these macros in selected locations. The
initial introduction of these macros simply replaces all identifiable nested 2D- and
3D-loops with single line statements (and adjusts indenting accordingly). Do loops
are identifiable if they comform to either:

DO jk = ....

DO jj = .... DO jj = ...

DO ji = .... DO ji = ...
. OR .
. .

END DO END DO

END DO END DO

END DO

and white-space variants thereof.

Additionally, only loops with recognised jj and ji loops limits are treated; these are:
Lower limits of 1, 2 or fs_2
Upper limits of jpi, jpim1 or fs_jpim1 (for ji) or jpj, jpjm1 or fs_jpjm1 (for jj)

The macro naming convention takes the form: DO_2D_BT_LR where:

B is the Bottom offset from the PE's inner domain;
T is the Top offset from the PE's inner domain;
L is the Left offset from the PE's inner domain;
R is the Right offset from the PE's inner domain

So, given an inner domain of 2,jpim1 and 2,jpjm1, a typical example would replace:

DO jj = 2, jpj

DO ji = 1, jpim1
.
.

END DO

END DO

with:

DO_2D_01_10
.
.
END_2D

similar conventions apply to the 3D loops macros. jk loop limits are retained
through macro arguments and are not restricted. This includes the possibility of
strides for which an extra set of DO_3DS macros are defined.

In the example definition below the inner PE domain is defined by start indices of
(kIs, kJs) and end indices of (kIe, KJe)

#define DO_2D_00_00 DO jj = kJs, kJe ; DO ji = kIs, kIe
#define END_2D END DO ; END DO

TO DO:


Only conventional nested loops have been identified and replaced by this step. There are constructs such as:

DO jk = 2, jpkm1

z2d(:,:) = z2d(:,:) + e3w(:,:,jk,Kmm) * z3d(:,:,jk) * wmask(:,:,jk)

END DO

which may need to be considered.

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