source: trunk/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zche.F90 @ 7646

Last change on this file since 7646 was 7646, checked in by timgraham, 5 years ago

Merge of dev_merge_2016 into trunk. UPDATE TO ARCHFILES NEEDED for XIOS2.
LIM_SRC_s/limrhg.F90 to follow in next commit due to change of kind (I'm unable to do it in this commit).
Merged using the following steps:

1) svn merge —reintegrate svn+ssh://forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/trunk .
2) Resolve minor conflicts in sette.sh and namelist_cfg for ORCA2LIM3 (due to a change in trunk after branch was created)
3) svn commit
4) svn switch svn+ssh://forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/trunk
5) svn merge svn+ssh://forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/branches/2016/dev_merge_2016 .
6) At this stage I checked out a clean copy of the branch to compare against what is about to be committed to the trunk.
6) svn commit #Commit code to the trunk

In this commit I have also reverted a change to Fcheck_archfile.sh which was causing problems on the Paris machine.

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