source: NEMO/branches/2020/dev_r12563_ASINTER-06_ABL_improvement/src/TOP/PISCES/P4Z/p4zche.F90

Last change on this file was 12377, checked in by acc, 10 months ago

The big one. Merging all 2019 developments from the option 1 branch back onto the trunk.

This changeset reproduces 2019/dev_r11943_MERGE_2019 on the trunk using a 2-URL merge
onto a working copy of the trunk. I.e.:

svn merge —ignore-ancestry \

svn+ssh://acc@forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/NEMO/trunk \
svn+ssh://acc@forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/NEMO/branches/2019/dev_r11943_MERGE_2019 ./

The —ignore-ancestry flag avoids problems that may otherwise arise from the fact that
the merge history been trunk and branch may have been applied in a different order but
care has been taken before this step to ensure that all applicable fixes and updates
are present in the merge branch.

The trunk state just before this step has been branched to releases/release-4.0-HEAD
and that branch has been immediately tagged as releases/release-4.0.2. Any fixes
or additions in response to tickets on 4.0, 4.0.1 or 4.0.2 should be done on
releases/release-4.0-HEAD. From now on future 'point' releases (e.g. 4.0.2) will
remain unchanged with periodic releases as needs demand. Note release-4.0-HEAD is a
transitional naming convention. Future full releases, say 4.2, will have a release-4.2
branch which fulfills this role and the first point release (e.g. 4.2.0) will be made
immediately following the release branch creation.

2020 developments can be started from any trunk revision later than this one.

  • Property svn:keywords set to Id
File size: 36.2 KB
Line 
1MODULE 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   !!----------------------------------------------------------------------
139CONTAINS
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.00023656*ztkel    &
200            &       + 0.0047036e-4*ztkel**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.0408*ztkel - 0.0327957*ztkel**2 + 0.0000316528*ztkel**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-6*gdept(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.2*ztr - 12.641 + 1.525*zisqrt   &
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.53*zsqrt - 77.942*zsal        &
284               &      + 1.728*zsal15 - 0.0996*zsal*zsal)*ztr         &
285               &      + (148.0248 + 137.1942*zsqrt + 1.62142*zsal)   &
286               &      + (-24.4344 - 25.085*zsqrt - 0.2474*zsal)      & 
287               &      * zlogt + 0.053105*zsqrt*ztkel
288
289               ! DISSOCIATION COEFFICIENT FOR CARBONATE ACCORDING TO
290               ! MEHRBACH (1973) REFIT BY MILLERO (1995), seawater scale
291               zck1    = -1.0*(3633.86*ztr - 61.2172 + 9.6777*zlogt  &
292                  - 0.011555*zsal + 0.0001152*zsal*zsal)
293               zck2    = -1.0*(471.78*ztr + 25.9290 - 3.16967*zlogt      &
294                  - 0.01781*zsal + 0.0001122*zsal*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.752*ztr + 115.540 - 18.453*zlogt   &
302              &          + (-106.736*ztr + 0.69171) * zsqrt       &
303              &          + (-0.65643*ztr - 0.01844) * zsal
304
305              zck2p    = -8814.715*ztr + 172.1033 - 27.927*zlogt  &
306              &          + (-160.340*ztr + 1.3566)*zsqrt          &
307              &          + (0.37335*ztr - 0.05778)*zsal
308
309              zck3p    = -3070.75*ztr - 18.126                    &
310              &          + (17.27039*ztr + 2.81197) * zsqrt       &
311              &          + (-44.99486*ztr - 0.09984) * zsal 
312
313              ! CONSTANT FOR SILICATE, MILLERO (1995)
314              zcksi    = -8904.2*ztr  + 117.400 - 19.334*zlogt   &
315              &          + (-458.79*ztr + 3.5913) * zisqrt       &
316              &          + (188.74*ztr - 1.5998) * zis           &
317              &          + (-12.1652*ztr + 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.077993*ztkel + 2839.319*ztr + 71.595*LOG10( ztkel )   &
323                  &      + (-0.77712 + 0.00284263*ztkel + 178.34*ztr) * zsqrt  &
324                  &      - 0.07711*zsal + 0.0041249*zsal15
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.3073*zis + 5254.0/ztkel1)
437               fesol(ji,jj,jk,2) = 10**(2.517 - 0.8885*zis**0.5 + 0.2139 * zis - 1320.0/ztkel1 )
438               fesol(ji,jj,jk,3) = 10**(0.4511 - 0.3305*zis**0.5 - 1996.0/ztkel1 )
439               fesol(ji,jj,jk,4) = 10**(-0.2965 - 0.7881*zis**0.5 - 4086.0/ztkel1 )
440               fesol(ji,jj,jk,5) = 10**(4.4466 - 0.8505*zis**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 = za2*za2 - 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) + zh*ak13(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/(zdeqndh*d[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]_old*10**(-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]_old*10**(-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/(zdeqndh*zh_prev)) -1| ~ |zeqn/(zdeqndh*zh_prev)|
760
761         ! Alternatively:
762         ! |\Delta pH| = |zeqn/(zdeqndh*zh_prev*LOG(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   !!======================================================================
820END MODULE p4zche
Note: See TracBrowser for help on using the repository browser.