Context Navigation

p4zche.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: 36.2 KB

Line
1	MODULE p4zche
2	!!======================================================================
3	!! * MODULE p4zche *
4	!! TOP : PISCES Sea water chemistry computed following OCMIP protocol
5	!!======================================================================
6	!! History : OPA ! 1988 (E. Maier-Reimer) Original code
7	!! - ! 1998 (O. Aumont) addition
8	!! - ! 1999 (C. Le Quere) modification
9	!! NEMO 1.0 ! 2004 (O. Aumont) modification
10	!! - ! 2006 (R. Gangsto) modification
11	!! 2.0 ! 2007-12 (C. Ethe, G. Madec) F90
12	!! ! 2011-02 (J. Simeon, J.Orr ) update O2 solubility constants
13	!! 3.6 ! 2016-03 (O. Aumont) Change chemistry to MOCSY standards
14	!!----------------------------------------------------------------------
15	!! p4z_che : Sea water chemistry computed following OCMIP protocol
16	!!----------------------------------------------------------------------
17	USE oce_trc ! shared variables between ocean and passive tracers
18	USE trc ! passive tracers common variables
19	USE sms_pisces ! PISCES Source Minus Sink variables
20	USE lib_mpp ! MPP library
21	USE eosbn2, ONLY : neos
22
23	IMPLICIT NONE
24	PRIVATE
25
26	PUBLIC p4z_che !
27	PUBLIC p4z_che_alloc !
28	PUBLIC ahini_for_at !
29	PUBLIC solve_at_general !
30
31	REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: sio3eq ! chemistry of Si
32	REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: fekeq ! chemistry of Fe
33	REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: chemc ! Solubilities of O2 and CO2
34	REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: chemo2 ! Solubilities of O2 and CO2
35	REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) :: fesol ! solubility of Fe
36	REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: salinprac ! Practical salinity
37	REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: tempis ! In situ temperature
38
39	REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: akb3 !: ???
40	REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: akw3 !: ???
41	REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: akf3 !: ???
42	REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: aks3 !: ???
43	REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: ak1p3 !: ???
44	REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: ak2p3 !: ???
45	REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: ak3p3 !: ???
46	REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: aksi3 !: ???
47	REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: borat !: ???
48	REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: fluorid !: ???
49	REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: sulfat !: ???
50
51	!!* Variable for chemistry of the CO2 cycle
52
53	REAL(wp), PUBLIC :: atcox = 0.20946 ! units atm
54
55	REAL(wp) :: o2atm = 1. / ( 1000. * 0.20946 )
56
57	REAL(wp) :: rgas = 83.14472 ! universal gas constants
58	REAL(wp) :: oxyco = 1. / 22.4144 ! converts from liters of an ideal gas to moles
59	! ! coeff. for seawater pressure correction : millero 95
60	! ! AGRIF doesn't like the DATA instruction
61	REAL(wp) :: devk10 = -25.5
62	REAL(wp) :: devk11 = -15.82
63	REAL(wp) :: devk12 = -29.48
64	REAL(wp) :: devk13 = -20.02
65	REAL(wp) :: devk14 = -18.03
66	REAL(wp) :: devk15 = -9.78
67	REAL(wp) :: devk16 = -48.76
68	REAL(wp) :: devk17 = -14.51
69	REAL(wp) :: devk18 = -23.12
70	REAL(wp) :: devk19 = -26.57
71	REAL(wp) :: devk110 = -29.48
72	!
73	REAL(wp) :: devk20 = 0.1271
74	REAL(wp) :: devk21 = -0.0219
75	REAL(wp) :: devk22 = 0.1622
76	REAL(wp) :: devk23 = 0.1119
77	REAL(wp) :: devk24 = 0.0466
78	REAL(wp) :: devk25 = -0.0090
79	REAL(wp) :: devk26 = 0.5304
80	REAL(wp) :: devk27 = 0.1211
81	REAL(wp) :: devk28 = 0.1758
82	REAL(wp) :: devk29 = 0.2020
83	REAL(wp) :: devk210 = 0.1622
84	!
85	REAL(wp) :: devk30 = 0.
86	REAL(wp) :: devk31 = 0.
87	REAL(wp) :: devk32 = 2.608E-3
88	REAL(wp) :: devk33 = -1.409e-3
89	REAL(wp) :: devk34 = 0.316e-3
90	REAL(wp) :: devk35 = -0.942e-3
91	REAL(wp) :: devk36 = 0.
92	REAL(wp) :: devk37 = -0.321e-3
93	REAL(wp) :: devk38 = -2.647e-3
94	REAL(wp) :: devk39 = -3.042e-3
95	REAL(wp) :: devk310 = -2.6080e-3
96	!
97	REAL(wp) :: devk40 = -3.08E-3
98	REAL(wp) :: devk41 = 1.13E-3
99	REAL(wp) :: devk42 = -2.84E-3
100	REAL(wp) :: devk43 = -5.13E-3
101	REAL(wp) :: devk44 = -4.53e-3
102	REAL(wp) :: devk45 = -3.91e-3
103	REAL(wp) :: devk46 = -11.76e-3
104	REAL(wp) :: devk47 = -2.67e-3
105	REAL(wp) :: devk48 = -5.15e-3
106	REAL(wp) :: devk49 = -4.08e-3
107	REAL(wp) :: devk410 = -2.84e-3
108	!
109	REAL(wp) :: devk50 = 0.0877E-3
110	REAL(wp) :: devk51 = -0.1475E-3
111	REAL(wp) :: devk52 = 0.
112	REAL(wp) :: devk53 = 0.0794E-3
113	REAL(wp) :: devk54 = 0.09e-3
114	REAL(wp) :: devk55 = 0.054e-3
115	REAL(wp) :: devk56 = 0.3692E-3
116	REAL(wp) :: devk57 = 0.0427e-3
117	REAL(wp) :: devk58 = 0.09e-3
118	REAL(wp) :: devk59 = 0.0714e-3
119	REAL(wp) :: devk510 = 0.0
120	!
121	! General parameters
122	REAL(wp), PARAMETER :: pp_rdel_ah_target = 1.E-4_wp
123	REAL(wp), PARAMETER :: pp_ln10 = 2.302585092994045684018_wp
124
125	! Maximum number of iterations for each method
126	INTEGER, PARAMETER :: jp_maxniter_atgen = 20
127
128	! Bookkeeping variables for each method
129	! - SOLVE_AT_GENERAL
130	INTEGER :: niter_atgen = jp_maxniter_atgen
131
132	!! * Substitutions
133	# include "do_loop_substitute.h90"
134	!!----------------------------------------------------------------------
135	!! NEMO/TOP 4.0 , NEMO Consortium (2018)
136	!! $Id$
137	!! Software governed by the CeCILL license (see ./LICENSE)
138	!!----------------------------------------------------------------------
139	CONTAINS
140
141	SUBROUTINE p4z_che( Kbb, Kmm )
142	!!---------------------------------------------------------------------
143	!! * ROUTINE p4z_che *
144	!!
145	!! ** Purpose : Sea water chemistry computed following OCMIP protocol
146	!!
147	!! ** Method : - ...
148	!!---------------------------------------------------------------------
149	INTEGER, INTENT(in) :: Kbb, Kmm ! time level indices
150	INTEGER :: ji, jj, jk
151	REAL(wp) :: ztkel, ztkel1, zt , zsal , zsal2 , zbuf1 , zbuf2
152	REAL(wp) :: ztgg , ztgg2, ztgg3 , ztgg4 , ztgg5
153	REAL(wp) :: zpres, ztc , zcl , zcpexp, zoxy , zcpexp2
154	REAL(wp) :: zsqrt, ztr , zlogt , zcek1, zc1, zplat
155	REAL(wp) :: zis , zis2 , zsal15, zisqrt, za1, za2
156	REAL(wp) :: zckb , zck1 , zck2 , zckw , zak1 , zak2 , zakb , zaksp0, zakw
157	REAL(wp) :: zck1p, zck2p, zck3p, zcksi, zak1p, zak2p, zak3p, zaksi
158	REAL(wp) :: zst , zft , zcks , zckf , zaksp1
159	REAL(wp) :: total2free, free2SWS, total2SWS, SWS2total
160	!!---------------------------------------------------------------------
161	!
162	IF( ln_timing ) CALL timing_start('p4z_che')
163	!
164	! Computation of chemical constants require practical salinity
165	! Thus, when TEOS08 is used, absolute salinity is converted to
166	! practical salinity
167	! -------------------------------------------------------------
168	IF (neos == -1) THEN
169	salinprac(:,:,:) = ts(:,:,:,jp_sal,Kmm) * 35.0 / 35.16504
170	ELSE
171	salinprac(:,:,:) = ts(:,:,:,jp_sal,Kmm)
172	ENDIF
173
174	!
175	! Computations of chemical constants require in situ temperature
176	! Here a quite simple formulation is used to convert
177	! potential temperature to in situ temperature. The errors is less than
178	! 0.04°C relative to an exact computation
179	! ---------------------------------------------------------------------
180	DO_3D_11_11( 1, jpk )
181	zpres = gdept(ji,jj,jk,Kmm) / 1000.
182	za1 = 0.04 * ( 1.0 + 0.185 * ts(ji,jj,jk,jp_tem,Kmm) + 0.035 * (salinprac(ji,jj,jk) - 35.0) )
183	za2 = 0.0075 * ( 1.0 - ts(ji,jj,jk,jp_tem,Kmm) / 30.0 )
184	tempis(ji,jj,jk) = ts(ji,jj,jk,jp_tem,Kmm) - za1 * zpres + za2 * zpres**2
185	END_3D
186	!
187	! CHEMICAL CONSTANTS - SURFACE LAYER
188	! ----------------------------------
189	!CDIR NOVERRCHK
190	DO jj = 1, jpj
191	!CDIR NOVERRCHK
192	DO ji = 1, jpi
193	! ! SET ABSOLUTE TEMPERATURE
194	ztkel = tempis(ji,jj,1) + 273.15
195	zt = ztkel * 0.01
196	zsal = salinprac(ji,jj,1) + ( 1.- tmask(ji,jj,1) ) * 35.
197	! ! LN(K0) OF SOLUBILITY OF CO2 (EQ. 12, WEISS, 1980)
198	! ! AND FOR THE ATMOSPHERE FOR NON IDEAL GAS
199	zcek1 = 9345.17/ztkel - 60.2409 + 23.3585 * LOG(zt) + zsal(0.023517 - 0.00023656ztkel &
200	& + 0.0047036e-4ztkel*2)
201	chemc(ji,jj,1) = EXP( zcek1 ) * 1E-6 * rhop(ji,jj,1) / 1000. ! mol/(L atm)
202	chemc(ji,jj,2) = -1636.75 + 12.0408ztkel - 0.0327957ztkel*2 + 0.0000316528ztkel**3
203	chemc(ji,jj,3) = 57.7 - 0.118*ztkel
204	!
205	END DO
206	END DO
207
208	! OXYGEN SOLUBILITY - DEEP OCEAN
209	! -------------------------------
210	!CDIR NOVERRCHK
211	DO jk = 1, jpk
212	!CDIR NOVERRCHK
213	DO jj = 1, jpj
214	!CDIR NOVERRCHK
215	DO ji = 1, jpi
216	ztkel = tempis(ji,jj,jk) + 273.15
217	zsal = salinprac(ji,jj,jk) + ( 1.- tmask(ji,jj,jk) ) * 35.
218	zsal2 = zsal * zsal
219	ztgg = LOG( ( 298.15 - tempis(ji,jj,jk) ) / ztkel ) ! Set the GORDON & GARCIA scaled temperature
220	ztgg2 = ztgg * ztgg
221	ztgg3 = ztgg2 * ztgg
222	ztgg4 = ztgg3 * ztgg
223	ztgg5 = ztgg4 * ztgg
224
225	zoxy = 2.00856 + 3.22400 * ztgg + 3.99063 * ztgg2 + 4.80299 * ztgg3 &
226	& + 9.78188e-1 * ztgg4 + 1.71069 * ztgg5 + zsal * ( -6.24097e-3 &
227	& - 6.93498e-3 * ztgg - 6.90358e-3 * ztgg2 - 4.29155e-3 * ztgg3 ) &
228	& - 3.11680e-7 * zsal2
229	chemo2(ji,jj,jk) = ( EXP( zoxy ) * o2atm ) * oxyco * atcox ! mol/(L atm)
230	END DO
231	END DO
232	END DO
233
234	! CHEMICAL CONSTANTS - DEEP OCEAN
235	! -------------------------------
236	!CDIR NOVERRCHK
237	DO jk = 1, jpk
238	!CDIR NOVERRCHK
239	DO jj = 1, jpj
240	!CDIR NOVERRCHK
241	DO ji = 1, jpi
242
243	! SET PRESSION ACCORDING TO SAUNDER (1980)
244	zplat = SIN ( ABS(gphit(ji,jj)*3.141592654/180.) )
245	zc1 = 5.92E-3 + zplat*2 5.25E-3
246	zpres = ((1-zc1)-SQRT(((1-zc1)*2)-(8.84E-6gdept(ji,jj,jk,Kmm)))) / 4.42E-6
247	zpres = zpres / 10.0
248
249	! SET ABSOLUTE TEMPERATURE
250	ztkel = tempis(ji,jj,jk) + 273.15
251	zsal = salinprac(ji,jj,jk) + ( 1.-tmask(ji,jj,jk) ) * 35.
252	zsqrt = SQRT( zsal )
253	zsal15 = zsqrt * zsal
254	zlogt = LOG( ztkel )
255	ztr = 1. / ztkel
256	zis = 19.924 * zsal / ( 1000.- 1.005 * zsal )
257	zis2 = zis * zis
258	zisqrt = SQRT( zis )
259	ztc = tempis(ji,jj,jk) + ( 1.- tmask(ji,jj,jk) ) * 20.
260
261	! CHLORINITY (WOOSTER ET AL., 1969)
262	zcl = zsal / 1.80655
263
264	! TOTAL SULFATE CONCENTR. [MOLES/kg soln]
265	zst = 0.14 * zcl /96.062
266
267	! TOTAL FLUORIDE CONCENTR. [MOLES/kg soln]
268	zft = 0.000067 * zcl /18.9984
269
270	! DISSOCIATION CONSTANT FOR SULFATES on free H scale (Dickson 1990)
271	zcks = EXP(-4276.1 * ztr + 141.328 - 23.093 * zlogt &
272	& + (-13856. * ztr + 324.57 - 47.986 * zlogt) * zisqrt &
273	& + (35474. * ztr - 771.54 + 114.723 * zlogt) * zis &
274	& - 2698. * ztr * zis*1.5 + 1776. ztr * zis2 &
275	& + LOG(1.0 - 0.001005 * zsal))
276
277	! DISSOCIATION CONSTANT FOR FLUORIDES on free H scale (Dickson and Riley 79)
278	zckf = EXP( 1590.2ztr - 12.641 + 1.525zisqrt &
279	& + LOG(1.0d0 - 0.001005d0*zsal) &
280	& + LOG(1.0d0 + zst/zcks))
281
282	! DISSOCIATION CONSTANT FOR CARBONATE AND BORATE
283	zckb= (-8966.90 - 2890.53zsqrt - 77.942zsal &
284	& + 1.728zsal15 - 0.0996zsalzsal)ztr &
285	& + (148.0248 + 137.1942zsqrt + 1.62142zsal) &
286	& + (-24.4344 - 25.085zsqrt - 0.2474zsal) &
287	& * zlogt + 0.053105zsqrtztkel
288
289	! DISSOCIATION COEFFICIENT FOR CARBONATE ACCORDING TO
290	! MEHRBACH (1973) REFIT BY MILLERO (1995), seawater scale
291	zck1 = -1.0(3633.86ztr - 61.2172 + 9.6777*zlogt &
292	- 0.011555zsal + 0.0001152zsal*zsal)
293	zck2 = -1.0(471.78ztr + 25.9290 - 3.16967*zlogt &
294	- 0.01781zsal + 0.0001122zsal*zsal)
295
296	! PKW (H2O) (MILLERO, 1995) from composite data
297	zckw = -13847.26 * ztr + 148.9652 - 23.6521 * zlogt + ( 118.67 * ztr &
298	- 5.977 + 1.0495 * zlogt ) * zsqrt - 0.01615 * zsal
299
300	! CONSTANTS FOR PHOSPHATE (MILLERO, 1995)
301	zck1p = -4576.752ztr + 115.540 - 18.453zlogt &
302	& + (-106.736ztr + 0.69171) zsqrt &
303	& + (-0.65643ztr - 0.01844) zsal
304
305	zck2p = -8814.715ztr + 172.1033 - 27.927zlogt &
306	& + (-160.340ztr + 1.3566)zsqrt &
307	& + (0.37335ztr - 0.05778)zsal
308
309	zck3p = -3070.75*ztr - 18.126 &
310	& + (17.27039ztr + 2.81197) zsqrt &
311	& + (-44.99486ztr - 0.09984) zsal
312
313	! CONSTANT FOR SILICATE, MILLERO (1995)
314	zcksi = -8904.2ztr + 117.400 - 19.334zlogt &
315	& + (-458.79ztr + 3.5913) zisqrt &
316	& + (188.74ztr - 1.5998) zis &
317	& + (-12.1652ztr + 0.07871) zis2 &
318	& + LOG(1.0 - 0.001005*zsal)
319
320	! APPARENT SOLUBILITY PRODUCT K'SP OF CALCITE IN SEAWATER
321	! (S=27-43, T=2-25 DEG C) at pres =0 (atmos. pressure) (MUCCI 1983)
322	zaksp0 = -171.9065 -0.077993ztkel + 2839.319ztr + 71.595*LOG10( ztkel ) &
323	& + (-0.77712 + 0.00284263ztkel + 178.34ztr) * zsqrt &
324	& - 0.07711zsal + 0.0041249zsal15
325
326	! CONVERT FROM DIFFERENT PH SCALES
327	total2free = 1.0/(1.0 + zst/zcks)
328	free2SWS = 1. + zst/zcks + zft/(zckf*total2free)
329	total2SWS = total2free * free2SWS
330	SWS2total = 1.0 / total2SWS
331
332	! K1, K2 OF CARBONIC ACID, KB OF BORIC ACID, KW (H2O) (LIT.?)
333	zak1 = 10*(zck1) total2SWS
334	zak2 = 10*(zck2) total2SWS
335	zakb = EXP( zckb ) * total2SWS
336	zakw = EXP( zckw )
337	zaksp1 = 10**(zaksp0)
338	zak1p = exp( zck1p )
339	zak2p = exp( zck2p )
340	zak3p = exp( zck3p )
341	zaksi = exp( zcksi )
342	zckf = zckf * total2SWS
343
344	! FORMULA FOR CPEXP AFTER EDMOND & GIESKES (1970)
345	! (REFERENCE TO CULBERSON & PYTKOQICZ (1968) AS MADE
346	! IN BROECKER ET AL. (1982) IS INCORRECT; HERE RGAS IS
347	! TAKEN TENFOLD TO CORRECT FOR THE NOTATION OF pres IN
348	! DBAR INSTEAD OF BAR AND THE EXPRESSION FOR CPEXP IS
349	! MULTIPLIED BY LN(10.) TO ALLOW USE OF EXP-FUNCTION
350	! WITH BASIS E IN THE FORMULA FOR AKSPP (CF. EDMOND
351	! & GIESKES (1970), P. 1285-1286 (THE SMALL
352	! FORMULA ON P. 1286 IS RIGHT AND CONSISTENT WITH THE
353	! SIGN IN PARTIAL MOLAR VOLUME CHANGE AS SHOWN ON P. 1285))
354	zcpexp = zpres / (rgas*ztkel)
355	zcpexp2 = zpres * zcpexp
356
357	! KB OF BORIC ACID, K1,K2 OF CARBONIC ACID PRESSURE
358	! CORRECTION AFTER CULBERSON AND PYTKOWICZ (1968)
359	! (CF. BROECKER ET AL., 1982)
360
361	zbuf1 = - ( devk10 + devk20 * ztc + devk30 * ztc * ztc )
362	zbuf2 = 0.5 * ( devk40 + devk50 * ztc )
363	ak13(ji,jj,jk) = zak1 * EXP( zbuf1 * zcpexp + zbuf2 * zcpexp2 )
364
365	zbuf1 = - ( devk11 + devk21 * ztc + devk31 * ztc * ztc )
366	zbuf2 = 0.5 * ( devk41 + devk51 * ztc )
367	ak23(ji,jj,jk) = zak2 * EXP( zbuf1 * zcpexp + zbuf2 * zcpexp2 )
368
369	zbuf1 = - ( devk12 + devk22 * ztc + devk32 * ztc * ztc )
370	zbuf2 = 0.5 * ( devk42 + devk52 * ztc )
371	akb3(ji,jj,jk) = zakb * EXP( zbuf1 * zcpexp + zbuf2 * zcpexp2 )
372
373	zbuf1 = - ( devk13 + devk23 * ztc + devk33 * ztc * ztc )
374	zbuf2 = 0.5 * ( devk43 + devk53 * ztc )
375	akw3(ji,jj,jk) = zakw * EXP( zbuf1 * zcpexp + zbuf2 * zcpexp2 )
376
377	zbuf1 = - ( devk14 + devk24 * ztc + devk34 * ztc * ztc )
378	zbuf2 = 0.5 * ( devk44 + devk54 * ztc )
379	aks3(ji,jj,jk) = zcks * EXP( zbuf1 * zcpexp + zbuf2 * zcpexp2 )
380
381	zbuf1 = - ( devk15 + devk25 * ztc + devk35 * ztc * ztc )
382	zbuf2 = 0.5 * ( devk45 + devk55 * ztc )
383	akf3(ji,jj,jk) = zckf * EXP( zbuf1 * zcpexp + zbuf2 * zcpexp2 )
384
385	zbuf1 = - ( devk17 + devk27 * ztc + devk37 * ztc * ztc )
386	zbuf2 = 0.5 * ( devk47 + devk57 * ztc )
387	ak1p3(ji,jj,jk) = zak1p * EXP( zbuf1 * zcpexp + zbuf2 * zcpexp2 )
388
389	zbuf1 = - ( devk18 + devk28 * ztc + devk38 * ztc * ztc )
390	zbuf2 = 0.5 * ( devk48 + devk58 * ztc )
391	ak2p3(ji,jj,jk) = zak2p * EXP( zbuf1 * zcpexp + zbuf2 * zcpexp2 )
392
393	zbuf1 = - ( devk19 + devk29 * ztc + devk39 * ztc * ztc )
394	zbuf2 = 0.5 * ( devk49 + devk59 * ztc )
395	ak3p3(ji,jj,jk) = zak3p * EXP( zbuf1 * zcpexp + zbuf2 * zcpexp2 )
396
397	zbuf1 = - ( devk110 + devk210 * ztc + devk310 * ztc * ztc )
398	zbuf2 = 0.5 * ( devk410 + devk510 * ztc )
399	aksi3(ji,jj,jk) = zaksi * EXP( zbuf1 * zcpexp + zbuf2 * zcpexp2 )
400
401	! CONVERT FROM DIFFERENT PH SCALES
402	total2free = 1.0/(1.0 + zst/aks3(ji,jj,jk))
403	free2SWS = 1. + zst/aks3(ji,jj,jk) + zft/akf3(ji,jj,jk)
404	total2SWS = total2free * free2SWS
405	SWS2total = 1.0 / total2SWS
406
407	! Convert to total scale
408	ak13(ji,jj,jk) = ak13(ji,jj,jk) * SWS2total
409	ak23(ji,jj,jk) = ak23(ji,jj,jk) * SWS2total
410	akb3(ji,jj,jk) = akb3(ji,jj,jk) * SWS2total
411	akw3(ji,jj,jk) = akw3(ji,jj,jk) * SWS2total
412	ak1p3(ji,jj,jk) = ak1p3(ji,jj,jk) * SWS2total
413	ak2p3(ji,jj,jk) = ak2p3(ji,jj,jk) * SWS2total
414	ak3p3(ji,jj,jk) = ak3p3(ji,jj,jk) * SWS2total
415	aksi3(ji,jj,jk) = aksi3(ji,jj,jk) * SWS2total
416	akf3(ji,jj,jk) = akf3(ji,jj,jk) / total2free
417
418	! APPARENT SOLUBILITY PRODUCT K'SP OF CALCITE
419	! AS FUNCTION OF PRESSURE FOLLOWING MILLERO
420	! (P. 1285) AND BERNER (1976)
421	zbuf1 = - ( devk16 + devk26 * ztc + devk36 * ztc * ztc )
422	zbuf2 = 0.5 * ( devk46 + devk56 * ztc )
423	aksp(ji,jj,jk) = zaksp1 * EXP( zbuf1 * zcpexp + zbuf2 * zcpexp2 )
424
425	! TOTAL F, S, and BORATE CONCENTR. [MOLES/L]
426	borat(ji,jj,jk) = 0.0002414 * zcl / 10.811
427	sulfat(ji,jj,jk) = zst
428	fluorid(ji,jj,jk) = zft
429
430	! Iron and SIO3 saturation concentration from ...
431	sio3eq(ji,jj,jk) = EXP( LOG( 10.) * ( 6.44 - 968. / ztkel ) ) * 1.e-6
432	fekeq (ji,jj,jk) = 10**( 17.27 - 1565.7 / ztkel )
433
434	! Liu and Millero (1999) only valid 5 - 50 degC
435	ztkel1 = MAX( 5. , tempis(ji,jj,jk) ) + 273.16
436	fesol(ji,jj,jk,1) = 10*(-13.486 - 0.1856 zis*0.5 + 0.3073zis + 5254.0/ztkel1)
437	fesol(ji,jj,jk,2) = 10*(2.517 - 0.8885zis*0.5 + 0.2139 zis - 1320.0/ztkel1 )
438	fesol(ji,jj,jk,3) = 10*(0.4511 - 0.3305zis**0.5 - 1996.0/ztkel1 )
439	fesol(ji,jj,jk,4) = 10*(-0.2965 - 0.7881zis**0.5 - 4086.0/ztkel1 )
440	fesol(ji,jj,jk,5) = 10*(4.4466 - 0.8505zis**0.5 - 7980.0/ztkel1 )
441	END DO
442	END DO
443	END DO
444	!
445	IF( ln_timing ) CALL timing_stop('p4z_che')
446	!
447	END SUBROUTINE p4z_che
448
449	SUBROUTINE ahini_for_at(p_hini, Kbb )
450	!!---------------------------------------------------------------------
451	!! * ROUTINE ahini_for_at *
452	!!
453	!! Subroutine returns the root for the 2nd order approximation of the
454	!! DIC -- B_T -- A_CB equation for [H+] (reformulated as a cubic
455	!! polynomial) around the local minimum, if it exists.
456	!! Returns * 1E-03_wp if p_alkcb <= 0
457	!! * 1E-10_wp if p_alkcb >= 2*p_dictot + p_bortot
458	!! * 1E-07_wp if 0 < p_alkcb < 2*p_dictot + p_bortot
459	!! and the 2nd order approximation does not have
460	!! a solution
461	!!---------------------------------------------------------------------
462	REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(OUT) :: p_hini
463	INTEGER, INTENT(in) :: Kbb ! time level indices
464	INTEGER :: ji, jj, jk
465	REAL(wp) :: zca1, zba1
466	REAL(wp) :: zd, zsqrtd, zhmin
467	REAL(wp) :: za2, za1, za0
468	REAL(wp) :: p_dictot, p_bortot, p_alkcb
469	!!---------------------------------------------------------------------
470
471	IF( ln_timing ) CALL timing_start('ahini_for_at')
472	!
473	DO_3D_11_11( 1, jpk )
474	p_alkcb = tr(ji,jj,jk,jptal,Kbb) * 1000. / (rhop(ji,jj,jk) + rtrn)
475	p_dictot = tr(ji,jj,jk,jpdic,Kbb) * 1000. / (rhop(ji,jj,jk) + rtrn)
476	p_bortot = borat(ji,jj,jk)
477	IF (p_alkcb <= 0.) THEN
478	p_hini(ji,jj,jk) = 1.e-3
479	ELSEIF (p_alkcb >= (2.*p_dictot + p_bortot)) THEN
480	p_hini(ji,jj,jk) = 1.e-10_wp
481	ELSE
482	zca1 = p_dictot/( p_alkcb + rtrn )
483	zba1 = p_bortot/ (p_alkcb + rtrn )
484	! Coefficients of the cubic polynomial
485	za2 = aKb3(ji,jj,jk)(1. - zba1) + ak13(ji,jj,jk)(1.-zca1)
486	za1 = ak13(ji,jj,jk)akb3(ji,jj,jk)(1. - zba1 - zca1) &
487	& + ak13(ji,jj,jk)ak23(ji,jj,jk)(1. - (zca1+zca1))
488	za0 = ak13(ji,jj,jk)ak23(ji,jj,jk)akb3(ji,jj,jk)*(1. - zba1 - (zca1+zca1))
489	! Taylor expansion around the minimum
490	zd = za2za2 - 3.za1 ! Discriminant of the quadratic equation
491	! for the minimum close to the root
492
493	IF(zd > 0.) THEN ! If the discriminant is positive
494	zsqrtd = SQRT(zd)
495	IF(za2 < 0) THEN
496	zhmin = (-za2 + zsqrtd)/3.
497	ELSE
498	zhmin = -za1/(za2 + zsqrtd)
499	ENDIF
500	p_hini(ji,jj,jk) = zhmin + SQRT(-(za0 + zhmin(za1 + zhmin(za2 + zhmin)))/zsqrtd)
501	ELSE
502	p_hini(ji,jj,jk) = 1.e-7
503	ENDIF
504	!
505	ENDIF
506	END_3D
507	!
508	IF( ln_timing ) CALL timing_stop('ahini_for_at')
509	!
510	END SUBROUTINE ahini_for_at
511
512	!===============================================================================
513
514	SUBROUTINE anw_infsup( p_alknw_inf, p_alknw_sup, Kbb )
515
516	! Subroutine returns the lower and upper bounds of "non-water-selfionization"
517	! contributions to total alkalinity (the infimum and the supremum), i.e
518	! inf(TA - [OH-] + [H+]) and sup(TA - [OH-] + [H+])
519
520	! Argument variables
521	REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(OUT) :: p_alknw_inf
522	REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(OUT) :: p_alknw_sup
523	INTEGER, INTENT(in) :: Kbb ! time level indices
524
525	p_alknw_inf(:,:,:) = -tr(:,:,:,jppo4,Kbb) * 1000. / (rhop(:,:,:) + rtrn) - sulfat(:,:,:) &
526	& - fluorid(:,:,:)
527	p_alknw_sup(:,:,:) = (2. * tr(:,:,:,jpdic,Kbb) + 2. * tr(:,:,:,jppo4,Kbb) + tr(:,:,:,jpsil,Kbb) ) &
528	& * 1000. / (rhop(:,:,:) + rtrn) + borat(:,:,:)
529
530	END SUBROUTINE anw_infsup
531
532
533	SUBROUTINE solve_at_general( p_hini, zhi, Kbb )
534
535	! Universal pH solver that converges from any given initial value,
536	! determines upper an lower bounds for the solution if required
537
538	! Argument variables
539	!--------------------
540	REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(IN) :: p_hini
541	REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(OUT) :: zhi
542	INTEGER, INTENT(in) :: Kbb ! time level indices
543
544	! Local variables
545	!-----------------
546	INTEGER :: ji, jj, jk, jn
547	REAL(wp) :: zh_ini, zh, zh_prev, zh_lnfactor
548	REAL(wp) :: zdelta, zh_delta
549	REAL(wp) :: zeqn, zdeqndh, zalka
550	REAL(wp) :: aphscale
551	REAL(wp) :: znumer_dic, zdnumer_dic, zdenom_dic, zalk_dic, zdalk_dic
552	REAL(wp) :: znumer_bor, zdnumer_bor, zdenom_bor, zalk_bor, zdalk_bor
553	REAL(wp) :: znumer_po4, zdnumer_po4, zdenom_po4, zalk_po4, zdalk_po4
554	REAL(wp) :: znumer_sil, zdnumer_sil, zdenom_sil, zalk_sil, zdalk_sil
555	REAL(wp) :: znumer_so4, zdnumer_so4, zdenom_so4, zalk_so4, zdalk_so4
556	REAL(wp) :: znumer_flu, zdnumer_flu, zdenom_flu, zalk_flu, zdalk_flu
557	REAL(wp) :: zalk_wat, zdalk_wat
558	REAL(wp) :: zfact, p_alktot, zdic, zbot, zpt, zst, zft, zsit
559	LOGICAL :: l_exitnow
560	REAL(wp), PARAMETER :: pz_exp_threshold = 1.0
561	REAL(wp), DIMENSION(jpi,jpj,jpk) :: zalknw_inf, zalknw_sup, rmask, zh_min, zh_max, zeqn_absmin
562
563	IF( ln_timing ) CALL timing_start('solve_at_general')
564
565	CALL anw_infsup( zalknw_inf, zalknw_sup, Kbb )
566
567	rmask(:,:,:) = tmask(:,:,:)
568	zhi(:,:,:) = 0.
569
570	! TOTAL H+ scale: conversion factor for Htot = aphscale * Hfree
571	DO_3D_11_11( 1, jpk )
572	IF (rmask(ji,jj,jk) == 1.) THEN
573	p_alktot = tr(ji,jj,jk,jptal,Kbb) * 1000. / (rhop(ji,jj,jk) + rtrn)
574	aphscale = 1. + sulfat(ji,jj,jk)/aks3(ji,jj,jk)
575	zh_ini = p_hini(ji,jj,jk)
576
577	zdelta = (p_alktot-zalknw_inf(ji,jj,jk))*2 + 4.akw3(ji,jj,jk)/aphscale
578
579	IF(p_alktot >= zalknw_inf(ji,jj,jk)) THEN
580	zh_min(ji,jj,jk) = 2.*akw3(ji,jj,jk) /( p_alktot-zalknw_inf(ji,jj,jk) + SQRT(zdelta) )
581	ELSE
582	zh_min(ji,jj,jk) = aphscale*(-(p_alktot-zalknw_inf(ji,jj,jk)) + SQRT(zdelta) ) / 2.
583	ENDIF
584
585	zdelta = (p_alktot-zalknw_sup(ji,jj,jk))*2 + 4.akw3(ji,jj,jk)/aphscale
586
587	IF(p_alktot <= zalknw_sup(ji,jj,jk)) THEN
588	zh_max(ji,jj,jk) = aphscale*(-(p_alktot-zalknw_sup(ji,jj,jk)) + SQRT(zdelta) ) / 2.
589	ELSE
590	zh_max(ji,jj,jk) = 2.*akw3(ji,jj,jk) /( p_alktot-zalknw_sup(ji,jj,jk) + SQRT(zdelta) )
591	ENDIF
592
593	zhi(ji,jj,jk) = MAX(MIN(zh_max(ji,jj,jk), zh_ini), zh_min(ji,jj,jk))
594	ENDIF
595	END_3D
596
597	zeqn_absmin(:,:,:) = HUGE(1._wp)
598
599	DO jn = 1, jp_maxniter_atgen
600	DO_3D_11_11( 1, jpk )
601	IF (rmask(ji,jj,jk) == 1.) THEN
602	zfact = rhop(ji,jj,jk) / 1000. + rtrn
603	p_alktot = tr(ji,jj,jk,jptal,Kbb) / zfact
604	zdic = tr(ji,jj,jk,jpdic,Kbb) / zfact
605	zbot = borat(ji,jj,jk)
606	zpt = tr(ji,jj,jk,jppo4,Kbb) / zfact * po4r
607	zsit = tr(ji,jj,jk,jpsil,Kbb) / zfact
608	zst = sulfat (ji,jj,jk)
609	zft = fluorid(ji,jj,jk)
610	aphscale = 1. + sulfat(ji,jj,jk)/aks3(ji,jj,jk)
611	zh = zhi(ji,jj,jk)
612	zh_prev = zh
613
614	! H2CO3 - HCO3 - CO3 : n=2, m=0
615	znumer_dic = 2.ak13(ji,jj,jk)ak23(ji,jj,jk) + zh*ak13(ji,jj,jk)
616	zdenom_dic = ak13(ji,jj,jk)ak23(ji,jj,jk) + zh(ak13(ji,jj,jk) + zh)
617	zalk_dic = zdic * (znumer_dic/zdenom_dic)
618	zdnumer_dic = ak13(ji,jj,jk)ak13(ji,jj,jk)ak23(ji,jj,jk) + zh &
619	(4.ak13(ji,jj,jk)ak23(ji,jj,jk) + zhak13(ji,jj,jk))
620	zdalk_dic = -zdic(zdnumer_dic/zdenom_dic*2)
621
622
623	! B(OH)3 - B(OH)4 : n=1, m=0
624	znumer_bor = akb3(ji,jj,jk)
625	zdenom_bor = akb3(ji,jj,jk) + zh
626	zalk_bor = zbot * (znumer_bor/zdenom_bor)
627	zdnumer_bor = akb3(ji,jj,jk)
628	zdalk_bor = -zbot(zdnumer_bor/zdenom_bor*2)
629
630
631	! H3PO4 - H2PO4 - HPO4 - PO4 : n=3, m=1
632	znumer_po4 = 3.ak1p3(ji,jj,jk)ak2p3(ji,jj,jk)*ak3p3(ji,jj,jk) &
633	& + zh(2.ak1p3(ji,jj,jk)ak2p3(ji,jj,jk) + zh ak1p3(ji,jj,jk))
634	zdenom_po4 = ak1p3(ji,jj,jk)ak2p3(ji,jj,jk)ak3p3(ji,jj,jk) &
635	& + zh( ak1p3(ji,jj,jk)ak2p3(ji,jj,jk) + zh*(ak1p3(ji,jj,jk) + zh))
636	zalk_po4 = zpt * (znumer_po4/zdenom_po4 - 1.) ! Zero level of H3PO4 = 1
637	zdnumer_po4 = ak1p3(ji,jj,jk)ak2p3(ji,jj,jk)ak1p3(ji,jj,jk)ak2p3(ji,jj,jk)ak3p3(ji,jj,jk) &
638	& + zh(4.ak1p3(ji,jj,jk)ak1p3(ji,jj,jk)ak2p3(ji,jj,jk)*ak3p3(ji,jj,jk) &
639	& + zh(9.ak1p3(ji,jj,jk)ak2p3(ji,jj,jk)ak3p3(ji,jj,jk) &
640	& + ak1p3(ji,jj,jk)ak1p3(ji,jj,jk)ak2p3(ji,jj,jk) &
641	& + zh(4.ak1p3(ji,jj,jk)ak2p3(ji,jj,jk) + zh ak1p3(ji,jj,jk) ) ) )
642	zdalk_po4 = -zpt * (zdnumer_po4/zdenom_po4**2)
643
644	! H4SiO4 - H3SiO4 : n=1, m=0
645	znumer_sil = aksi3(ji,jj,jk)
646	zdenom_sil = aksi3(ji,jj,jk) + zh
647	zalk_sil = zsit * (znumer_sil/zdenom_sil)
648	zdnumer_sil = aksi3(ji,jj,jk)
649	zdalk_sil = -zsit * (zdnumer_sil/zdenom_sil**2)
650
651	! HSO4 - SO4 : n=1, m=1
652	aphscale = 1.0 + zst/aks3(ji,jj,jk)
653	znumer_so4 = aks3(ji,jj,jk) * aphscale
654	zdenom_so4 = aks3(ji,jj,jk) * aphscale + zh
655	zalk_so4 = zst * (znumer_so4/zdenom_so4 - 1.)
656	zdnumer_so4 = aks3(ji,jj,jk)
657	zdalk_so4 = -zst * (zdnumer_so4/zdenom_so4**2)
658
659	! HF - F : n=1, m=1
660	znumer_flu = akf3(ji,jj,jk)
661	zdenom_flu = akf3(ji,jj,jk) + zh
662	zalk_flu = zft * (znumer_flu/zdenom_flu - 1.)
663	zdnumer_flu = akf3(ji,jj,jk)
664	zdalk_flu = -zft * (zdnumer_flu/zdenom_flu**2)
665
666	! H2O - OH
667	aphscale = 1.0 + zst/aks3(ji,jj,jk)
668	zalk_wat = akw3(ji,jj,jk)/zh - zh/aphscale
669	zdalk_wat = -akw3(ji,jj,jk)/zh**2 - 1./aphscale
670
671	! CALCULATE [ALK]([CO3--], [HCO3-])
672	zeqn = zalk_dic + zalk_bor + zalk_po4 + zalk_sil &
673	& + zalk_so4 + zalk_flu &
674	& + zalk_wat - p_alktot
675
676	zalka = p_alktot - (zalk_bor + zalk_po4 + zalk_sil &
677	& + zalk_so4 + zalk_flu + zalk_wat)
678
679	zdeqndh = zdalk_dic + zdalk_bor + zdalk_po4 + zdalk_sil &
680	& + zdalk_so4 + zdalk_flu + zdalk_wat
681
682	! Adapt bracketing interval
683	IF(zeqn > 0._wp) THEN
684	zh_min(ji,jj,jk) = zh_prev
685	ELSEIF(zeqn < 0._wp) THEN
686	zh_max(ji,jj,jk) = zh_prev
687	ENDIF
688
689	IF(ABS(zeqn) >= 0.5_wp*zeqn_absmin(ji,jj,jk)) THEN
690	! if the function evaluation at the current point is
691	! not decreasing faster than with a bisection step (at least linearly)
692	! in absolute value take one bisection step on [ph_min, ph_max]
693	! ph_new = (ph_min + ph_max)/2d0
694	!
695	! In terms of [H]_new:
696	! [H]_new = 10**(-ph_new)
697	! = 10**(-(ph_min + ph_max)/2d0)
698	! = SQRT(10**(-(ph_min + phmax)))
699	! = SQRT(zh_max * zh_min)
700	zh = SQRT(zh_max(ji,jj,jk) * zh_min(ji,jj,jk))
701	zh_lnfactor = (zh - zh_prev)/zh_prev ! Required to test convergence below
702	ELSE
703	! dzeqn/dpH = dzeqn/d[H] * d[H]/dpH
704	! = -zdeqndh * LOG(10) * [H]
705	! \Delta pH = -zeqn/(zdeqndhd[H]/dpH) = zeqn/(zdeqndh[H]*LOG(10))
706	!
707	! pH_new = pH_old + \deltapH
708	!
709	! [H]_new = 10**(-pH_new)
710	! = 10**(-pH_old - \Delta pH)
711	! = [H]_old * 10*(-zeqn/(zdeqndh[H]_old*LOG(10)))
712	! = [H]_old * EXP(-LOG(10)zeqn/(zdeqndh[H]_old*LOG(10)))
713	! = [H]_old * EXP(-zeqn/(zdeqndh*[H]_old))
714
715	zh_lnfactor = -zeqn/(zdeqndh*zh_prev)
716
717	IF(ABS(zh_lnfactor) > pz_exp_threshold) THEN
718	zh = zh_prev*EXP(zh_lnfactor)
719	ELSE
720	zh_delta = zh_lnfactor*zh_prev
721	zh = zh_prev + zh_delta
722	ENDIF
723
724	IF( zh < zh_min(ji,jj,jk) ) THEN
725	! if [H]_new < [H]_min
726	! i.e., if ph_new > ph_max then
727	! take one bisection step on [ph_prev, ph_max]
728	! ph_new = (ph_prev + ph_max)/2d0
729	! In terms of [H]_new:
730	! [H]_new = 10**(-ph_new)
731	! = 10**(-(ph_prev + ph_max)/2d0)
732	! = SQRT(10**(-(ph_prev + phmax)))
733	! = SQRT([H]_old10*(-ph_max))
734	! = SQRT([H]_old * zh_min)
735	zh = SQRT(zh_prev * zh_min(ji,jj,jk))
736	zh_lnfactor = (zh - zh_prev)/zh_prev ! Required to test convergence below
737	ENDIF
738
739	IF( zh > zh_max(ji,jj,jk) ) THEN
740	! if [H]_new > [H]_max
741	! i.e., if ph_new < ph_min, then
742	! take one bisection step on [ph_min, ph_prev]
743	! ph_new = (ph_prev + ph_min)/2d0
744	! In terms of [H]_new:
745	! [H]_new = 10**(-ph_new)
746	! = 10**(-(ph_prev + ph_min)/2d0)
747	! = SQRT(10**(-(ph_prev + ph_min)))
748	! = SQRT([H]_old10*(-ph_min))
749	! = SQRT([H]_old * zhmax)
750	zh = SQRT(zh_prev * zh_max(ji,jj,jk))
751	zh_lnfactor = (zh - zh_prev)/zh_prev ! Required to test convergence below
752	ENDIF
753	ENDIF
754
755	zeqn_absmin(ji,jj,jk) = MIN( ABS(zeqn), zeqn_absmin(ji,jj,jk))
756
757	! Stop iterations once \|\delta{[H]}/[H]\| < rdel
758	! <=> \|(zh - zh_prev)/zh_prev\| = \|EXP(-zeqn/(zdeqndh*zh_prev)) -1\| < rdel
759	! \|EXP(-zeqn/(zdeqndhzh_prev)) -1\| ~ \|zeqn/(zdeqndhzh_prev)\|
760
761	! Alternatively:
762	! \|\Delta pH\| = \|zeqn/(zdeqndhzh_prevLOG(10))\|
763	! ~ 1/LOG(10) * \|\Delta [H]\|/[H]
764	! < 1/LOG(10) * rdel
765
766	! Hence \|zeqn/(zdeqndh*zh)\| < rdel
767
768	! rdel <-- pp_rdel_ah_target
769	l_exitnow = (ABS(zh_lnfactor) < pp_rdel_ah_target)
770
771	IF(l_exitnow) THEN
772	rmask(ji,jj,jk) = 0.
773	ENDIF
774
775	zhi(ji,jj,jk) = zh
776
777	IF(jn >= jp_maxniter_atgen) THEN
778	zhi(ji,jj,jk) = -1._wp
779	ENDIF
780
781	ENDIF
782	END_3D
783	END DO
784	!
785
786	IF( ln_timing ) CALL timing_stop('solve_at_general')
787	!
788	END SUBROUTINE solve_at_general
789
790
791	INTEGER FUNCTION p4z_che_alloc()
792	!!----------------------------------------------------------------------
793	!! * ROUTINE p4z_che_alloc *
794	!!----------------------------------------------------------------------
795	INTEGER :: ierr(3) ! Local variables
796	!!----------------------------------------------------------------------
797
798	ierr(:) = 0
799
800	ALLOCATE( sio3eq(jpi,jpj,jpk), fekeq(jpi,jpj,jpk), chemc(jpi,jpj,3), chemo2(jpi,jpj,jpk), STAT=ierr(1) )
801
802	ALLOCATE( akb3(jpi,jpj,jpk) , tempis(jpi, jpj, jpk), &
803	& akw3(jpi,jpj,jpk) , borat (jpi,jpj,jpk) , &
804	& aks3(jpi,jpj,jpk) , akf3(jpi,jpj,jpk) , &
805	& ak1p3(jpi,jpj,jpk) , ak2p3(jpi,jpj,jpk) , &
806	& ak3p3(jpi,jpj,jpk) , aksi3(jpi,jpj,jpk) , &
807	& fluorid(jpi,jpj,jpk) , sulfat(jpi,jpj,jpk) , &
808	& salinprac(jpi,jpj,jpk), STAT=ierr(2) )
809
810	ALLOCATE( fesol(jpi,jpj,jpk,5), STAT=ierr(3) )
811
812	!* Variable for chemistry of the CO2 cycle
813	p4z_che_alloc = MAXVAL( ierr )
814	!
815	IF( p4z_che_alloc /= 0 ) CALL ctl_stop( 'STOP', 'p4z_che_alloc : failed to allocate arrays.' )
816	!
817	END FUNCTION p4z_che_alloc
818
819	!!======================================================================
820	END MODULE p4zche

Note: See TracBrowser for help on using the repository browser.

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

Context Navigation

source: NEMO/branches/2019/dev_r11943_MERGE_2019/src/TOP/PISCES/P4Z/p4zche.F90 @ 12340

Download in other formats: