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.

p4zche.F90 in branches/2015/nemo_v3_6_STABLE/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z – NEMO

Context Navigation

source: branches/2015/nemo_v3_6_STABLE/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zche.F90 @ 7607

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

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

Download in other formats: