Index: branches/2015/nemo_v3_6_STABLE/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zche.F90
===================================================================
--- branches/2015/nemo_v3_6_STABLE/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zche.F90	(revision 7522)
+++ branches/2015/nemo_v3_6_STABLE/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zche.F90	(revision 7607)
@@ -11,8 +11,9 @@
    !!             2.0  !  2007-12  (C. Ethe, G. Madec)  F90
    !!                  !  2011-02  (J. Simeon, J.Orr ) update O2 solubility constants
+   !!             3.6  !  2016-03  (O. Aumont) Change chemistry to MOCSY standards
    !!----------------------------------------------------------------------
 #if defined key_pisces
    !!----------------------------------------------------------------------
-   !!   'key_pisces'                                       PISCES bio-model
+   !!   'key_pisces*'                                      PISCES bio-model
    !!----------------------------------------------------------------------
    !!   p4z_che      :  Sea water chemistry computed following OCMIP protocol
@@ -22,78 +23,115 @@
    USE sms_pisces    !  PISCES Source Minus Sink variables
    USE lib_mpp       !  MPP library
+   USE eosbn2, ONLY : nn_eos
 
    IMPLICIT NONE
    PRIVATE
 
-   PUBLIC   p4z_che         !
-   PUBLIC   p4z_che_alloc   !
-
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   sio3eq   ! chemistry of Si
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   fekeq    ! chemistry of Fe
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   chemc    ! Solubilities of O2 and CO2
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   chemo2   ! Solubilities of O2 and CO2
+   PUBLIC   p4z_che          !
+   PUBLIC   p4z_che_alloc    !
+   PUBLIC   p4z_che_ahini    !
+   PUBLIC   p4z_che_solve_hi !
+
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)   :: sio3eq   ! chemistry of Si
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)   :: fekeq    ! chemistry of Fe
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)   :: chemc    ! Solubilities of O2 and CO2
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)   :: chemo2    ! Solubilities of O2 and CO2
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) :: fesol    ! solubility of Fe
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   tempis   ! In situ temperature
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   salinprac  ! Practical salinity
+
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   akb3       !: ???
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   akw3       !: ???
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   akf3       !: ???
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   aks3       !: ???
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   ak1p3      !: ???
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   ak2p3      !: ???
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   ak3p3      !: ???
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   aksi3      !: ???
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   borat      !: ???
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   fluorid    !: ???
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   sulfat     !: ???
+
+   !!* Variable for chemistry of the CO2 cycle
 
    REAL(wp), PUBLIC ::   atcox  = 0.20946         ! units atm
 
-   REAL(wp) ::   salchl = 1. / 1.80655    ! conversion factor for salinity --> chlorinity (Wooster et al. 1969)
    REAL(wp) ::   o2atm  = 1. / ( 1000. * 0.20946 )  
 
-   REAL(wp) ::   rgas   = 83.14472       ! universal gas constants
-   REAL(wp) ::   oxyco  = 1. / 22.4144   ! converts from liters of an ideal gas to moles
-
-   REAL(wp) ::   bor1   = 0.00023        ! borat constants
-   REAL(wp) ::   bor2   = 1. / 10.82
-
-   REAL(wp) ::   st1    =      0.14     ! constants for calculate concentrations for sulfate
-   REAL(wp) ::   st2    =  1./96.062    !  (Morris & Riley 1966)
-
-   REAL(wp) ::   ft1    =    0.000067   ! constants for calculate concentrations for fluorides
-   REAL(wp) ::   ft2    = 1./18.9984    ! (Dickson & Riley 1979 )
-
-   !                                    ! volumetric solubility constants for o2 in ml/L  
-   REAL(wp) ::   ox0    =  2.00856      ! from Table 1 for Eq 8 of Garcia and Gordon, 1992.
-   REAL(wp) ::   ox1    =  3.22400      ! corrects for moisture and fugacity, but not total atmospheric pressure
-   REAL(wp) ::   ox2    =  3.99063      !      Original PISCES code noted this was a solubility, but 
-   REAL(wp) ::   ox3    =  4.80299      ! was in fact a bunsen coefficient with units L-O2/(Lsw atm-O2)
-   REAL(wp) ::   ox4    =  9.78188e-1   ! Hence, need to divide EXP( zoxy ) by 1000, ml-O2 => L-O2
-   REAL(wp) ::   ox5    =  1.71069      ! and atcox = 0.20946 to add the 1/atm dimension.
-   REAL(wp) ::   ox6    = -6.24097e-3   
-   REAL(wp) ::   ox7    = -6.93498e-3 
-   REAL(wp) ::   ox8    = -6.90358e-3
-   REAL(wp) ::   ox9    = -4.29155e-3 
-   REAL(wp) ::   ox10   = -3.11680e-7 
+   REAL(wp) ::   rgas   = 83.14472      ! universal gas constants
+   REAL(wp) ::   oxyco  = 1. / 22.4144  ! converts from liters of an ideal gas to moles
 
    !                                    ! coeff. for seawater pressure correction : millero 95
    !                                    ! AGRIF doesn't like the DATA instruction
-   REAL(wp) :: devk11  = -25.5
-   REAL(wp) :: devk12  = -15.82
-   REAL(wp) :: devk13  = -29.48
-   REAL(wp) :: devk14  = -25.60
-   REAL(wp) :: devk15  = -48.76
+   REAL(wp) :: devk10  = -25.5
+   REAL(wp) :: devk11  = -15.82
+   REAL(wp) :: devk12  = -29.48
+   REAL(wp) :: devk13  = -20.02
+   REAL(wp) :: devk14  = -18.03
+   REAL(wp) :: devk15  = -9.78
+   REAL(wp) :: devk16  = -48.76
+   REAL(wp) :: devk17  = -14.51
+   REAL(wp) :: devk18  = -23.12
+   REAL(wp) :: devk19  = -26.57
+   REAL(wp) :: devk110  = -29.48
    !
-   REAL(wp) :: devk21  = 0.1271
-   REAL(wp) :: devk22  = -0.0219
-   REAL(wp) :: devk23  = 0.1622
-   REAL(wp) :: devk24  = 0.2324
-   REAL(wp) :: devk25  = 0.5304
+   REAL(wp) :: devk20  = 0.1271
+   REAL(wp) :: devk21  = -0.0219
+   REAL(wp) :: devk22  = 0.1622
+   REAL(wp) :: devk23  = 0.1119
+   REAL(wp) :: devk24  = 0.0466
+   REAL(wp) :: devk25  = -0.0090
+   REAL(wp) :: devk26  = 0.5304
+   REAL(wp) :: devk27  = 0.1211
+   REAL(wp) :: devk28  = 0.1758
+   REAL(wp) :: devk29  = 0.2020
+   REAL(wp) :: devk210  = 0.1622
    !
+   REAL(wp) :: devk30  = 0.
    REAL(wp) :: devk31  = 0.
-   REAL(wp) :: devk32  = 0.
-   REAL(wp) :: devk33  = 2.608E-3
-   REAL(wp) :: devk34  = -3.6246E-3
-   REAL(wp) :: devk35  = 0.
+   REAL(wp) :: devk32  = 2.608E-3
+   REAL(wp) :: devk33  = -1.409e-3
+   REAL(wp) :: devk34  = 0.316e-3
+   REAL(wp) :: devk35  = -0.942e-3
+   REAL(wp) :: devk36  = 0.
+   REAL(wp) :: devk37  = -0.321e-3
+   REAL(wp) :: devk38  = -2.647e-3
+   REAL(wp) :: devk39  = -3.042e-3
+   REAL(wp) :: devk310  = -2.6080e-3
    !
-   REAL(wp) :: devk41  = -3.08E-3
-   REAL(wp) :: devk42  = 1.13E-3
-   REAL(wp) :: devk43  = -2.84E-3
-   REAL(wp) :: devk44  = -5.13E-3
-   REAL(wp) :: devk45  = -11.76E-3
+   REAL(wp) :: devk40  = -3.08E-3
+   REAL(wp) :: devk41  = 1.13E-3
+   REAL(wp) :: devk42  = -2.84E-3
+   REAL(wp) :: devk43  = -5.13E-3
+   REAL(wp) :: devk44  = -4.53e-3
+   REAL(wp) :: devk45  = -3.91e-3
+   REAL(wp) :: devk46  = -11.76e-3
+   REAL(wp) :: devk47  = -2.67e-3
+   REAL(wp) :: devk48  = -5.15e-3
+   REAL(wp) :: devk49  = -4.08e-3
+   REAL(wp) :: devk410  = -2.84e-3
    !
-   REAL(wp) :: devk51  = 0.0877E-3
-   REAL(wp) :: devk52  = -0.1475E-3     
-   REAL(wp) :: devk53  = 0.
-   REAL(wp) :: devk54  = 0.0794E-3      
-   REAL(wp) :: devk55  = 0.3692E-3      
+   REAL(wp) :: devk50  = 0.0877E-3
+   REAL(wp) :: devk51  = -0.1475E-3     
+   REAL(wp) :: devk52  = 0.
+   REAL(wp) :: devk53  = 0.0794E-3      
+   REAL(wp) :: devk54  = 0.09e-3
+   REAL(wp) :: devk55  = 0.054e-3
+   REAL(wp) :: devk56  = 0.3692E-3
+   REAL(wp) :: devk57  = 0.0427e-3
+   REAL(wp) :: devk58  = 0.09e-3
+   REAL(wp) :: devk59  = 0.0714e-3
+   REAL(wp) :: devk510  = 0.0
+   !
+   ! General parameters
+   REAL(wp), PARAMETER :: pp_rdel_ah_target = 1.E-4_wp
+   REAL(wp), PARAMETER :: pp_ln10 = 2.302585092994045684018_wp
+
+   ! Maximum number of iterations for each method
+   INTEGER, PARAMETER :: jp_maxniter_atgen    = 20
+
+   ! Bookkeeping variables for each method
+   ! - SOLVE_AT_GENERAL
+   INTEGER :: niter_atgen    = jp_maxniter_atgen
 
    !!* Substitution
@@ -115,14 +153,28 @@
       !!---------------------------------------------------------------------
       INTEGER  ::   ji, jj, jk
-      REAL(wp) ::   ztkel, zt   , zt2   , zsal  , zsal2 , zbuf1 , zbuf2
+      REAL(wp) ::   ztkel, ztkel1, zt , zsal  , zsal2 , zbuf1 , zbuf2
       REAL(wp) ::   ztgg , ztgg2, ztgg3 , ztgg4 , ztgg5
       REAL(wp) ::   zpres, ztc  , zcl   , zcpexp, zoxy  , zcpexp2
       REAL(wp) ::   zsqrt, ztr  , zlogt , zcek1, zc1, zplat
-      REAL(wp) ::   zis  , zis2 , zsal15, zisqrt, za1  , za2
+      REAL(wp) ::   zis  , zis2 , zsal15, zisqrt, za1, za2
       REAL(wp) ::   zckb , zck1 , zck2  , zckw  , zak1 , zak2  , zakb , zaksp0, zakw
+      REAL(wp) ::   zck1p, zck2p, zck3p, zcksi, zak1p, zak2p, zak3p, zaksi
       REAL(wp) ::   zst  , zft  , zcks  , zckf  , zaksp1
+      REAL(wp) ::   total2free, free2SWS, total2SWS, SWS2total
+
       !!---------------------------------------------------------------------
       !
       IF( nn_timing == 1 )  CALL timing_start('p4z_che')
+      !
+      ! Computation of chemical constants require practical salinity
+      ! Thus, when TEOS08 is used, absolute salinity is converted to 
+      ! practical salinity
+      ! -------------------------------------------------------------
+      IF (nn_eos == -1) THEN
+         salinprac(:,:,:) = tsn(:,:,:,jp_sal) * 35.0 / 35.16504
+      ELSE
+         salinprac(:,:,:) = tsn(:,:,:,jp_sal)
+      ENDIF
+
       !
       ! Computations of chemical constants require in situ temperature
@@ -135,5 +187,5 @@
             DO ji = 1, jpi
                zpres = fsdept(ji,jj,jk) / 1000.
-               za1 = 0.04 * ( 1.0 + 0.185 * tsn(ji,jj,jk,jp_tem) + 0.035 * (tsn(ji,jj,jk,jp_sal) - 35.0) )
+               za1 = 0.04 * ( 1.0 + 0.185 * tsn(ji,jj,jk,jp_tem) + 0.035 * (salinprac(ji,jj,jk) - 35.0) )
                za2 = 0.0075 * ( 1.0 - tsn(ji,jj,jk,jp_tem) / 30.0 )
                tempis(ji,jj,jk) = tsn(ji,jj,jk,jp_tem) - za1 * zpres + za2 * zpres**2
@@ -151,14 +203,10 @@
             ztkel = tempis(ji,jj,1) + 273.15
             zt    = ztkel * 0.01
-            zt2   = zt * zt
-            zsal  = tsn(ji,jj,1,jp_sal) + ( 1.- tmask(ji,jj,1) ) * 35.
-            zsal2 = zsal * zsal
-            zlogt = LOG( zt )
+            zsal  = salinprac(ji,jj,1) + ( 1.- tmask(ji,jj,1) ) * 35.
             !                             ! LN(K0) OF SOLUBILITY OF CO2 (EQ. 12, WEISS, 1980)
             !                             !     AND FOR THE ATMOSPHERE FOR NON IDEAL GAS
             zcek1 = 9345.17/ztkel - 60.2409 + 23.3585 * LOG(zt) + zsal*(0.023517 - 0.00023656*ztkel    &
             &       + 0.0047036e-4*ztkel**2)
-            !                             ! SET SOLUBILITIES OF O2 AND CO2 
-            chemc(ji,jj,1) = EXP( zcek1 ) * 1.e-6 * rhop(ji,jj,1) / 1000. ! mol/(kg uatm)
+            chemc(ji,jj,1) = EXP( zcek1 ) * 1E-6 * rhop(ji,jj,1) / 1000. ! mol/(L atm)
             chemc(ji,jj,2) = -1636.75 + 12.0408*ztkel - 0.0327957*ztkel**2 + 0.0000316528*ztkel**3
             chemc(ji,jj,3) = 57.7 - 0.118*ztkel
@@ -176,5 +224,5 @@
             DO ji = 1, jpi
               ztkel = tempis(ji,jj,jk) + 273.15
-              zsal  = tsn(ji,jj,jk,jp_sal) + ( 1.- tmask(ji,jj,jk) ) * 35.
+              zsal  = salinprac(ji,jj,jk) + ( 1.- tmask(ji,jj,jk) ) * 35.
               zsal2 = zsal * zsal
               ztgg  = LOG( ( 298.15 - tempis(ji,jj,jk) ) / ztkel )  ! Set the GORDON & GARCIA scaled temperature
@@ -183,6 +231,9 @@
               ztgg4 = ztgg3 * ztgg
               ztgg5 = ztgg4 * ztgg
-              zoxy  = ox0 + ox1 * ztgg + ox2 * ztgg2 + ox3 * ztgg3 + ox4 * ztgg4 + ox5 * ztgg5   &
-                     + zsal * ( ox6 + ox7 * ztgg + ox8 * ztgg2 + ox9 * ztgg3 ) +  ox10 * zsal2
+
+              zoxy  = 2.00856 + 3.22400 * ztgg + 3.99063 * ztgg2 + 4.80299 * ztgg3    &
+              &       + 9.78188e-1 * ztgg4 + 1.71069 * ztgg5 + zsal * ( -6.24097e-3   &
+              &       - 6.93498e-3 * ztgg - 6.90358e-3 * ztgg2 - 4.29155e-3 * ztgg3 )   &
+              &       - 3.11680e-7 * zsal2
               chemo2(ji,jj,jk) = ( EXP( zoxy ) * o2atm ) * oxyco * atcox     ! mol/(L atm)
             END DO
@@ -209,5 +260,5 @@
                ! SET ABSOLUTE TEMPERATURE
                ztkel   = tempis(ji,jj,jk) + 273.15
-               zsal    = tsn(ji,jj,jk,jp_sal) + ( 1.-tmask(ji,jj,jk) ) * 35.
+               zsal    = salinprac(ji,jj,jk) + ( 1.-tmask(ji,jj,jk) ) * 35.
                zsqrt  = SQRT( zsal )
                zsal15  = zsqrt * zsal
@@ -220,11 +271,11 @@
 
                ! CHLORINITY (WOOSTER ET AL., 1969)
-               zcl     = zsal * salchl
+               zcl     = zsal / 1.80655
 
                ! TOTAL SULFATE CONCENTR. [MOLES/kg soln]
-               zst     = st1 * zcl * st2
+               zst     = 0.14 * zcl /96.062
 
                ! TOTAL FLUORIDE CONCENTR. [MOLES/kg soln]
-               zft     = ft1 * zcl * ft2
+               zft     = 0.000067 * zcl /18.9984
 
                ! DISSOCIATION CONSTANT FOR SULFATES on free H scale (Dickson 1990)
@@ -234,6 +285,4 @@
                &         - 2698. * ztr * zis**1.5 + 1776.* ztr * zis2         &
                &         + LOG(1.0 - 0.001005 * zsal))
-               !
-               aphscale(ji,jj,jk) = ( 1. + zst / zcks )
 
                ! DISSOCIATION CONSTANT FOR FLUORIDES on free H scale (Dickson and Riley 79)
@@ -249,5 +298,4 @@
                &      * zlogt + 0.053105*zsqrt*ztkel
 
-
                ! DISSOCIATION COEFFICIENT FOR CARBONATE ACCORDING TO 
                ! MEHRBACH (1973) REFIT BY MILLERO (1995), seawater scale
@@ -257,8 +305,27 @@
                   - 0.01781*zsal + 0.0001122*zsal*zsal)
 
-               ! PKW (H2O) (DICKSON AND RILEY, 1979)
-               zckw = -13847.26*ztr + 148.9652 - 23.6521 * zlogt    & 
-               &     + (118.67*ztr - 5.977 + 1.0495 * zlogt)        &
-               &     * zsqrt - 0.01615 * zsal
+               ! PKW (H2O) (MILLERO, 1995) from composite data
+               zckw    = -13847.26 * ztr + 148.9652 - 23.6521 * zlogt + ( 118.67 * ztr    &
+                         - 5.977 + 1.0495 * zlogt ) * zsqrt - 0.01615 * zsal
+
+               ! CONSTANTS FOR PHOSPHATE (MILLERO, 1995)
+              zck1p    = -4576.752*ztr + 115.540 - 18.453*zlogt   &
+              &          + (-106.736*ztr + 0.69171) * zsqrt       &
+              &          + (-0.65643*ztr - 0.01844) * zsal
+
+              zck2p    = -8814.715*ztr + 172.1033 - 27.927*zlogt  &
+              &          + (-160.340*ztr + 1.3566)*zsqrt          &
+              &          + (0.37335*ztr - 0.05778)*zsal
+
+              zck3p    = -3070.75*ztr - 18.126                    &
+              &          + (17.27039*ztr + 2.81197) * zsqrt       &
+              &          + (-44.99486*ztr - 0.09984) * zsal 
+
+              ! CONSTANT FOR SILICATE, MILLERO (1995)
+              zcksi    = -8904.2*ztr  + 117.400 - 19.334*zlogt   &
+              &          + (-458.79*ztr + 3.5913) * zisqrt       &
+              &          + (188.74*ztr - 1.5998) * zis           &
+              &          + (-12.1652*ztr + 0.07871) * zis2       &
+              &          + LOG(1.0 - 0.001005*zsal)
 
                ! APPARENT SOLUBILITY PRODUCT K'SP OF CALCITE IN SEAWATER
@@ -268,10 +335,21 @@
                   &      - 0.07711*zsal + 0.0041249*zsal15
 
+               ! CONVERT FROM DIFFERENT PH SCALES
+               total2free  = 1.0/(1.0 + zst/zcks)
+               free2SWS    = 1. + zst/zcks + zft/(zckf*total2free)
+               total2SWS   = total2free * free2SWS
+               SWS2total   = 1.0 / total2SWS
+
                ! K1, K2 OF CARBONIC ACID, KB OF BORIC ACID, KW (H2O) (LIT.?)
-               zak1    = 10**(zck1)
-               zak2    = 10**(zck2)
-               zakb    = EXP( zckb  )
+               zak1    = 10**(zck1) * total2SWS
+               zak2    = 10**(zck2) * total2SWS
+               zakb    = EXP( zckb ) * total2SWS
                zakw    = EXP( zckw )
                zaksp1  = 10**(zaksp0)
+               zak1p   = exp( zck1p )
+               zak2p   = exp( zck2p )
+               zak3p   = exp( zck3p )
+               zaksi   = exp( zcksi )
+               zckf    = zckf * total2SWS
 
                ! FORMULA FOR CPEXP AFTER EDMOND & GIESKES (1970)
@@ -285,6 +363,6 @@
                !        FORMULA ON P. 1286 IS RIGHT AND CONSISTENT WITH THE
                !        SIGN IN PARTIAL MOLAR VOLUME CHANGE AS SHOWN ON P. 1285))
-               zcpexp  = zpres /(rgas*ztkel)
-               zcpexp2 = zpres * zpres/(rgas*ztkel)
+               zcpexp  = zpres / (rgas*ztkel)
+               zcpexp2 = zpres * zcpexp
 
                ! KB OF BORIC ACID, K1,K2 OF CARBONIC ACID PRESSURE
@@ -292,35 +370,84 @@
                !        (CF. BROECKER ET AL., 1982)
 
-               zbuf1  = -     ( devk11 + devk21 * ztc + devk31 * ztc * ztc )
+               zbuf1  = -     ( devk10 + devk20 * ztc + devk30 * ztc * ztc )
+               zbuf2  = 0.5 * ( devk40 + devk50 * ztc )
+               ak13(ji,jj,jk) = zak1 * EXP( zbuf1 * zcpexp + zbuf2 * zcpexp2 )
+
+               zbuf1  =     - ( devk11 + devk21 * ztc + devk31 * ztc * ztc )
                zbuf2  = 0.5 * ( devk41 + devk51 * ztc )
-               ak13(ji,jj,jk) = zak1 * EXP( zbuf1 * zcpexp + zbuf2 * zcpexp2 )
+               ak23(ji,jj,jk) = zak2 * EXP( zbuf1 * zcpexp + zbuf2 * zcpexp2 )
 
                zbuf1  =     - ( devk12 + devk22 * ztc + devk32 * ztc * ztc )
                zbuf2  = 0.5 * ( devk42 + devk52 * ztc )
-               ak23(ji,jj,jk) = zak2 * EXP( zbuf1 * zcpexp + zbuf2 * zcpexp2 )
+               akb3(ji,jj,jk) = zakb * EXP( zbuf1 * zcpexp + zbuf2 * zcpexp2 )
 
                zbuf1  =     - ( devk13 + devk23 * ztc + devk33 * ztc * ztc )
                zbuf2  = 0.5 * ( devk43 + devk53 * ztc )
-               akb3(ji,jj,jk) = zakb * EXP( zbuf1 * zcpexp + zbuf2 * zcpexp2 )
+               akw3(ji,jj,jk) = zakw * EXP( zbuf1 * zcpexp + zbuf2 * zcpexp2 )
 
                zbuf1  =     - ( devk14 + devk24 * ztc + devk34 * ztc * ztc )
                zbuf2  = 0.5 * ( devk44 + devk54 * ztc )
-               akw3(ji,jj,jk) = zakw * EXP( zbuf1 * zcpexp + zbuf2 * zcpexp2 )
-
+               aks3(ji,jj,jk) = zcks * EXP( zbuf1 * zcpexp + zbuf2 * zcpexp2 )
+
+               zbuf1  =     - ( devk15 + devk25 * ztc + devk35 * ztc * ztc )
+               zbuf2  = 0.5 * ( devk45 + devk55 * ztc )
+               akf3(ji,jj,jk) = zckf * EXP( zbuf1 * zcpexp + zbuf2 * zcpexp2 )
+
+               zbuf1  =     - ( devk17 + devk27 * ztc + devk37 * ztc * ztc )
+               zbuf2  = 0.5 * ( devk47 + devk57 * ztc )
+               ak1p3(ji,jj,jk) = zak1p * EXP( zbuf1 * zcpexp + zbuf2 * zcpexp2 )
+
+               zbuf1  =     - ( devk18 + devk28 * ztc + devk38 * ztc * ztc )
+               zbuf2  = 0.5 * ( devk48 + devk58 * ztc )
+               ak2p3(ji,jj,jk) = zak2p * EXP( zbuf1 * zcpexp + zbuf2 * zcpexp2 )
+
+               zbuf1  =     - ( devk19 + devk29 * ztc + devk39 * ztc * ztc )
+               zbuf2  = 0.5 * ( devk49 + devk59 * ztc )
+               ak3p3(ji,jj,jk) = zak3p * EXP( zbuf1 * zcpexp + zbuf2 * zcpexp2 )
+
+               zbuf1  =     - ( devk110 + devk210 * ztc + devk310 * ztc * ztc )
+               zbuf2  = 0.5 * ( devk410 + devk510 * ztc )
+               aksi3(ji,jj,jk) = zaksi * EXP( zbuf1 * zcpexp + zbuf2 * zcpexp2 )
+
+               ! CONVERT FROM DIFFERENT PH SCALES
+               total2free  = 1.0/(1.0 + zst/aks3(ji,jj,jk))
+               free2SWS    = 1. + zst/aks3(ji,jj,jk) + zft/akf3(ji,jj,jk)
+               total2SWS   = total2free * free2SWS
+               SWS2total   = 1.0 / total2SWS
+
+               ! Convert to total scale
+               ak13(ji,jj,jk)  = ak13(ji,jj,jk)  * SWS2total
+               ak23(ji,jj,jk)  = ak23(ji,jj,jk)  * SWS2total
+               akb3(ji,jj,jk)  = akb3(ji,jj,jk)  * SWS2total
+               akw3(ji,jj,jk)  = akw3(ji,jj,jk)  * SWS2total
+               ak1p3(ji,jj,jk) = ak1p3(ji,jj,jk) * SWS2total
+               ak2p3(ji,jj,jk) = ak2p3(ji,jj,jk) * SWS2total
+               ak3p3(ji,jj,jk) = ak3p3(ji,jj,jk) * SWS2total
+               aksi3(ji,jj,jk) = aksi3(ji,jj,jk) * SWS2total
+               akf3(ji,jj,jk)  = akf3(ji,jj,jk)  / total2free
 
                ! APPARENT SOLUBILITY PRODUCT K'SP OF CALCITE 
                !        AS FUNCTION OF PRESSURE FOLLOWING MILLERO
                !        (P. 1285) AND BERNER (1976)
-               zbuf1  =     - ( devk15 + devk25 * ztc + devk35 * ztc * ztc )
-               zbuf2  = 0.5 * ( devk45 + devk55 * ztc )
+               zbuf1  =     - ( devk16 + devk26 * ztc + devk36 * ztc * ztc )
+               zbuf2  = 0.5 * ( devk46 + devk56 * ztc )
                aksp(ji,jj,jk) = zaksp1 * EXP( zbuf1 * zcpexp + zbuf2 * zcpexp2 )
 
-               ! TOTAL BORATE CONCENTR. [MOLES/L]
-               borat(ji,jj,jk) = bor1 * zcl * bor2
+               ! TOTAL F, S, and BORATE CONCENTR. [MOLES/L]
+               borat(ji,jj,jk) = 0.0002414 * zcl / 10.811
+               sulfat(ji,jj,jk) = zst
+               fluorid(ji,jj,jk) = zft 
 
                ! Iron and SIO3 saturation concentration from ...
                sio3eq(ji,jj,jk) = EXP(  LOG( 10.) * ( 6.44 - 968. / ztkel )  ) * 1.e-6
-               fekeq (ji,jj,jk) = 10**( 17.27 - 1565.7 / ( 273.15 + ztc ) )
-
+               fekeq (ji,jj,jk) = 10**( 17.27 - 1565.7 / ztkel )
+
+               ! Liu and Millero (1999) only valid 5 - 50 degC
+               ztkel1 = MAX( 5. , tempis(ji,jj,jk) ) + 273.16
+               fesol(ji,jj,jk,1) = 10**((-13.486) - (0.1856* (zis**0.5)) + (0.3073*zis) + (5254.0/ztkel1))
+               fesol(ji,jj,jk,2) = 10**(2.517 - (0.885*(zis**0.5)) + (0.2139 * zis) - (1320.0/ztkel1) )
+               fesol(ji,jj,jk,3) = 10**(0.4511 - (0.3305*(zis**0.5)) - (1996.0/ztkel1) )
+               fesol(ji,jj,jk,4) = 10**(-0.2965 - (0.7881*(zis**0.5)) - (4086.0/ztkel1) )
+               fesol(ji,jj,jk,5) = 10**(4.4466 - (0.8505*(zis**0.5)) - (7980.0/ztkel1) )
             END DO
          END DO
@@ -331,4 +458,358 @@
    END SUBROUTINE p4z_che
 
+   SUBROUTINE p4z_che_ahini( p_hini )
+      !!---------------------------------------------------------------------
+      !!                     ***  ROUTINE ahini_for_at  ***
+      !!
+      !! Subroutine returns the root for the 2nd order approximation of the
+      !! DIC -- B_T -- A_CB equation for [H+] (reformulated as a cubic 
+      !! polynomial) around the local minimum, if it exists.
+      !! Returns * 1E-03_wp if p_alkcb <= 0
+      !!         * 1E-10_wp if p_alkcb >= 2*p_dictot + p_bortot
+      !!         * 1E-07_wp if 0 < p_alkcb < 2*p_dictot + p_bortot
+      !!                    and the 2nd order approximation does not have 
+      !!                    a solution
+      !!---------------------------------------------------------------------
+      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(OUT)  ::  p_hini
+      INTEGER  ::   ji, jj, jk
+      REAL(wp)  ::  zca1, zba1
+      REAL(wp)  ::  zd, zsqrtd, zhmin
+      REAL(wp)  ::  za2, za1, za0
+      REAL(wp)  ::  p_dictot, p_bortot, p_alkcb 
+
+      IF( nn_timing == 1 )  CALL timing_start('p4z_che_ahini')
+      !
+      DO jk = 1, jpk
+        DO jj = 1, jpj
+          DO ji = 1, jpi
+            p_alkcb  = trb(ji,jj,jk,jptal) * 1000. / (rhop(ji,jj,jk) + rtrn)
+            p_dictot = trb(ji,jj,jk,jpdic) * 1000. / (rhop(ji,jj,jk) + rtrn)
+            p_bortot = borat(ji,jj,jk)
+            IF (p_alkcb <= 0.) THEN
+                p_hini(ji,jj,jk) = 1.e-3
+            ELSEIF (p_alkcb >= (2.*p_dictot + p_bortot)) THEN
+                p_hini(ji,jj,jk) = 1.e-10_wp
+            ELSE
+                zca1 = p_dictot/( p_alkcb + rtrn )
+                zba1 = p_bortot/ (p_alkcb + rtrn )
+           ! Coefficients of the cubic polynomial
+                za2 = aKb3(ji,jj,jk)*(1. - zba1) + ak13(ji,jj,jk)*(1.-zca1)
+                za1 = ak13(ji,jj,jk)*akb3(ji,jj,jk)*(1. - zba1 - zca1)    &
+                &     + ak13(ji,jj,jk)*ak23(ji,jj,jk)*(1. - (zca1+zca1))
+                za0 = ak13(ji,jj,jk)*ak23(ji,jj,jk)*akb3(ji,jj,jk)*(1. - zba1 - (zca1+zca1))
+                                        ! Taylor expansion around the minimum
+                zd = za2*za2 - 3.*za1   ! Discriminant of the quadratic equation
+                                        ! for the minimum close to the root
+
+                IF(zd > 0.) THEN        ! If the discriminant is positive
+                  zsqrtd = SQRT(zd)
+                  IF(za2 < 0) THEN
+                    zhmin = (-za2 + zsqrtd)/3.
+                  ELSE
+                    zhmin = -za1/(za2 + zsqrtd)
+                  ENDIF
+                  p_hini(ji,jj,jk) = zhmin + SQRT(-(za0 + zhmin*(za1 + zhmin*(za2 + zhmin)))/zsqrtd)
+                ELSE
+                  p_hini(ji,jj,jk) = 1.e-7
+                ENDIF
+             !
+             ENDIF
+          END DO
+        END DO
+      END DO
+      !
+      IF( nn_timing == 1 )  CALL timing_stop('p4z_che_ahini')
+      !
+   END SUBROUTINE p4z_che_ahini
+
+   !===============================================================================
+   SUBROUTINE anw_infsup( p_alknw_inf, p_alknw_sup )
+
+   ! Subroutine returns the lower and upper bounds of "non-water-selfionization"
+   ! contributions to total alkalinity (the infimum and the supremum), i.e
+   ! inf(TA - [OH-] + [H+]) and sup(TA - [OH-] + [H+])
+
+   ! Argument variables
+   REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(OUT) :: p_alknw_inf
+   REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(OUT) :: p_alknw_sup
+
+   p_alknw_inf(:,:,:) =  -trb(:,:,:,jppo4) * 1000. / (rhop(:,:,:) + rtrn) - sulfat(:,:,:)  &
+   &              - fluorid(:,:,:)
+   p_alknw_sup(:,:,:) =   (2. * trb(:,:,:,jpdic) + 2. * trb(:,:,:,jppo4) + trb(:,:,:,jpsil) )    &
+   &               * 1000. / (rhop(:,:,:) + rtrn) + borat(:,:,:) 
+
+   END SUBROUTINE anw_infsup
+
+
+   SUBROUTINE p4z_che_solve_hi( p_hini, zhi )
+
+   ! Universal pH solver that converges from any given initial value,
+   ! determines upper an lower bounds for the solution if required
+
+   ! Argument variables
+   !--------------------
+   REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(IN)   :: p_hini
+   REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(OUT)  :: zhi
+
+   ! Local variables
+   !-----------------
+   INTEGER   ::  ji, jj, jk, jn
+   REAL(wp)  ::  zh_ini, zh, zh_prev, zh_lnfactor
+   REAL(wp)  ::  zdelta, zh_delta
+   REAL(wp)  ::  zeqn, zdeqndh, zalka
+   REAL(wp)  ::  aphscale
+   REAL(wp)  ::  znumer_dic, zdnumer_dic, zdenom_dic, zalk_dic, zdalk_dic
+   REAL(wp)  ::  znumer_bor, zdnumer_bor, zdenom_bor, zalk_bor, zdalk_bor
+   REAL(wp)  ::  znumer_po4, zdnumer_po4, zdenom_po4, zalk_po4, zdalk_po4
+   REAL(wp)  ::  znumer_sil, zdnumer_sil, zdenom_sil, zalk_sil, zdalk_sil
+   REAL(wp)  ::  znumer_so4, zdnumer_so4, zdenom_so4, zalk_so4, zdalk_so4
+   REAL(wp)  ::  znumer_flu, zdnumer_flu, zdenom_flu, zalk_flu, zdalk_flu
+   REAL(wp)  ::  zalk_wat, zdalk_wat
+   REAL(wp)  ::  zfact, p_alktot, zdic, zbot, zpt, zst, zft, zsit
+   LOGICAL   ::  l_exitnow
+   REAL(wp), PARAMETER :: pz_exp_threshold = 1.0
+   REAL(wp), POINTER, DIMENSION(:,:,:) :: zalknw_inf, zalknw_sup, rmask, zh_min, zh_max, zeqn_absmin
+
+   IF( nn_timing == 1 )  CALL timing_start('p4z_che_solve_hi')
+      !  Allocate temporary workspace
+   CALL wrk_alloc( jpi, jpj, jpk, zalknw_inf, zalknw_sup, rmask )
+   CALL wrk_alloc( jpi, jpj, jpk, zh_min, zh_max, zeqn_absmin )
+
+   CALL anw_infsup( zalknw_inf, zalknw_sup )
+
+   rmask(:,:,:) = tmask(:,:,:)
+   zhi(:,:,:)   = 0.
+
+   ! TOTAL H+ scale: conversion factor for Htot = aphscale * Hfree
+   DO jk = 1, jpk
+      DO jj = 1, jpj
+         DO ji = 1, jpi
+            IF (rmask(ji,jj,jk) == 1.) THEN
+               p_alktot = trb(ji,jj,jk,jptal) * 1000. / (rhop(ji,jj,jk) + rtrn)
+               aphscale = 1. + sulfat(ji,jj,jk)/aks3(ji,jj,jk)
+               zh_ini = p_hini(ji,jj,jk)
+
+               zdelta = (p_alktot-zalknw_inf(ji,jj,jk))**2 + 4.*akw3(ji,jj,jk)/aphscale
+
+               IF(p_alktot >= zalknw_inf(ji,jj,jk)) THEN
+                 zh_min(ji,jj,jk) = 2.*akw3(ji,jj,jk) /( p_alktot-zalknw_inf(ji,jj,jk) + SQRT(zdelta) )
+               ELSE
+                 zh_min(ji,jj,jk) = aphscale*(-(p_alktot-zalknw_inf(ji,jj,jk)) + SQRT(zdelta) ) / 2.
+               ENDIF
+
+               zdelta = (p_alktot-zalknw_sup(ji,jj,jk))**2 + 4.*akw3(ji,jj,jk)/aphscale
+
+               IF(p_alktot <= zalknw_sup(ji,jj,jk)) THEN
+                 zh_max(ji,jj,jk) = aphscale*(-(p_alktot-zalknw_sup(ji,jj,jk)) + SQRT(zdelta) ) / 2.
+               ELSE
+                 zh_max(ji,jj,jk) = 2.*akw3(ji,jj,jk) /( p_alktot-zalknw_sup(ji,jj,jk) + SQRT(zdelta) )
+               ENDIF
+
+               zhi(ji,jj,jk) = MAX(MIN(zh_max(ji,jj,jk), zh_ini), zh_min(ji,jj,jk))
+            ENDIF
+         END DO
+      END DO
+   END DO
+
+   zeqn_absmin(:,:,:) = HUGE(1._wp)
+
+   DO jn = 1, jp_maxniter_atgen 
+   DO jk = 1, jpk
+      DO jj = 1, jpj
+         DO ji = 1, jpi
+            IF (rmask(ji,jj,jk) == 1.) THEN
+               zfact = rhop(ji,jj,jk) / 1000. + rtrn
+               p_alktot = trb(ji,jj,jk,jptal) / zfact
+               zdic  = trb(ji,jj,jk,jpdic) / zfact
+               zbot  = borat(ji,jj,jk)
+               zpt = trb(ji,jj,jk,jppo4) / zfact * po4r
+               zsit = trb(ji,jj,jk,jpsil) / zfact
+               zst = sulfat (ji,jj,jk)
+               zft = fluorid(ji,jj,jk)
+               aphscale = 1. + sulfat(ji,jj,jk)/aks3(ji,jj,jk)
+               zh = zhi(ji,jj,jk)
+               zh_prev = zh
+
+               ! H2CO3 - HCO3 - CO3 : n=2, m=0
+               znumer_dic = 2.*ak13(ji,jj,jk)*ak23(ji,jj,jk) + zh*ak13(ji,jj,jk)
+               zdenom_dic = ak13(ji,jj,jk)*ak23(ji,jj,jk) + zh*(ak13(ji,jj,jk) + zh)
+               zalk_dic   = zdic * (znumer_dic/zdenom_dic)
+               zdnumer_dic = ak13(ji,jj,jk)*ak13(ji,jj,jk)*ak23(ji,jj,jk) + zh     &
+                             *(4.*ak13(ji,jj,jk)*ak23(ji,jj,jk) + zh*ak13(ji,jj,jk))
+               zdalk_dic   = -zdic*(zdnumer_dic/zdenom_dic**2)
+
+
+               ! B(OH)3 - B(OH)4 : n=1, m=0
+               znumer_bor = akb3(ji,jj,jk)
+               zdenom_bor = akb3(ji,jj,jk) + zh
+               zalk_bor   = zbot * (znumer_bor/zdenom_bor)
+               zdnumer_bor = akb3(ji,jj,jk)
+               zdalk_bor   = -zbot*(zdnumer_bor/zdenom_bor**2)
+
+
+               ! H3PO4 - H2PO4 - HPO4 - PO4 : n=3, m=1
+               znumer_po4 = 3.*ak1p3(ji,jj,jk)*ak2p3(ji,jj,jk)*ak3p3(ji,jj,jk)  &
+               &            + zh*(2.*ak1p3(ji,jj,jk)*ak2p3(ji,jj,jk) + zh* ak1p3(ji,jj,jk))
+               zdenom_po4 = ak1p3(ji,jj,jk)*ak2p3(ji,jj,jk)*ak3p3(ji,jj,jk)     &
+               &            + zh*( ak1p3(ji,jj,jk)*ak2p3(ji,jj,jk) + zh*(ak1p3(ji,jj,jk) + zh))
+               zalk_po4   = zpt * (znumer_po4/zdenom_po4 - 1.) ! Zero level of H3PO4 = 1
+               zdnumer_po4 = ak1p3(ji,jj,jk)*ak2p3(ji,jj,jk)*ak1p3(ji,jj,jk)*ak2p3(ji,jj,jk)*ak3p3(ji,jj,jk)  &
+               &             + zh*(4.*ak1p3(ji,jj,jk)*ak1p3(ji,jj,jk)*ak2p3(ji,jj,jk)*ak3p3(ji,jj,jk)         &
+               &             + zh*(9.*ak1p3(ji,jj,jk)*ak2p3(ji,jj,jk)*ak3p3(ji,jj,jk)                         &
+               &             + ak1p3(ji,jj,jk)*ak1p3(ji,jj,jk)*ak2p3(ji,jj,jk)                                &
+               &             + zh*(4.*ak1p3(ji,jj,jk)*ak2p3(ji,jj,jk) + zh * ak1p3(ji,jj,jk) ) ) )
+               zdalk_po4   = -zpt * (zdnumer_po4/zdenom_po4**2)
+
+               ! H4SiO4 - H3SiO4 : n=1, m=0
+               znumer_sil = aksi3(ji,jj,jk)
+               zdenom_sil = aksi3(ji,jj,jk) + zh
+               zalk_sil   = zsit * (znumer_sil/zdenom_sil)
+               zdnumer_sil = aksi3(ji,jj,jk)
+               zdalk_sil   = -zsit * (zdnumer_sil/zdenom_sil**2)
+
+               ! HSO4 - SO4 : n=1, m=1
+               aphscale = 1.0 + zst/aks3(ji,jj,jk)
+               znumer_so4 = aks3(ji,jj,jk) * aphscale
+               zdenom_so4 = aks3(ji,jj,jk) * aphscale + zh
+               zalk_so4   = zst * (znumer_so4/zdenom_so4 - 1.)
+               zdnumer_so4 = aks3(ji,jj,jk)
+               zdalk_so4   = -zst * (zdnumer_so4/zdenom_so4**2)
+
+               ! HF - F : n=1, m=1
+               znumer_flu =  akf3(ji,jj,jk)
+               zdenom_flu =  akf3(ji,jj,jk) + zh
+               zalk_flu   =  zft * (znumer_flu/zdenom_flu - 1.)
+               zdnumer_flu = akf3(ji,jj,jk)
+               zdalk_flu   = -zft * (zdnumer_flu/zdenom_flu**2)
+
+               ! H2O - OH
+               aphscale = 1.0 + zst/aks3(ji,jj,jk)
+               zalk_wat   = akw3(ji,jj,jk)/zh - zh/aphscale
+               zdalk_wat  = -akw3(ji,jj,jk)/zh**2 - 1./aphscale
+
+               ! CALCULATE [ALK]([CO3--], [HCO3-])
+               zeqn = zalk_dic + zalk_bor + zalk_po4 + zalk_sil   &
+               &      + zalk_so4 + zalk_flu                       &
+               &      + zalk_wat - p_alktot
+
+               zalka = p_alktot - (zalk_bor + zalk_po4 + zalk_sil   &
+               &       + zalk_so4 + zalk_flu + zalk_wat)
+
+               zdeqndh = zdalk_dic + zdalk_bor + zdalk_po4 + zdalk_sil &
+               &         + zdalk_so4 + zdalk_flu + zdalk_wat
+
+               ! Adapt bracketing interval
+               IF(zeqn > 0._wp) THEN
+                 zh_min(ji,jj,jk) = zh_prev
+               ELSEIF(zeqn < 0._wp) THEN
+                 zh_max(ji,jj,jk) = zh_prev
+               ENDIF
+
+               IF(ABS(zeqn) >= 0.5_wp*zeqn_absmin(ji,jj,jk)) THEN
+               ! if the function evaluation at the current point is
+               ! not decreasing faster than with a bisection step (at least linearly)
+               ! in absolute value take one bisection step on [ph_min, ph_max]
+               ! ph_new = (ph_min + ph_max)/2d0
+               !
+               ! In terms of [H]_new:
+               ! [H]_new = 10**(-ph_new)
+               !         = 10**(-(ph_min + ph_max)/2d0)
+               !         = SQRT(10**(-(ph_min + phmax)))
+               !         = SQRT(zh_max * zh_min)
+                  zh = SQRT(zh_max(ji,jj,jk) * zh_min(ji,jj,jk))
+                  zh_lnfactor = (zh - zh_prev)/zh_prev ! Required to test convergence below
+               ELSE
+               ! dzeqn/dpH = dzeqn/d[H] * d[H]/dpH
+               !           = -zdeqndh * LOG(10) * [H]
+               ! \Delta pH = -zeqn/(zdeqndh*d[H]/dpH) = zeqn/(zdeqndh*[H]*LOG(10))
+               !
+               ! pH_new = pH_old + \deltapH
+               !
+               ! [H]_new = 10**(-pH_new)
+               !         = 10**(-pH_old - \Delta pH)
+               !         = [H]_old * 10**(-zeqn/(zdeqndh*[H]_old*LOG(10)))
+               !         = [H]_old * EXP(-LOG(10)*zeqn/(zdeqndh*[H]_old*LOG(10)))
+               !         = [H]_old * EXP(-zeqn/(zdeqndh*[H]_old))
+
+                  zh_lnfactor = -zeqn/(zdeqndh*zh_prev)
+
+                  IF(ABS(zh_lnfactor) > pz_exp_threshold) THEN
+                     zh          = zh_prev*EXP(zh_lnfactor)
+                  ELSE
+                     zh_delta    = zh_lnfactor*zh_prev
+                     zh          = zh_prev + zh_delta
+                  ENDIF
+
+                  IF( zh < zh_min(ji,jj,jk) ) THEN
+                     ! if [H]_new < [H]_min
+                     ! i.e., if ph_new > ph_max then
+                     ! take one bisection step on [ph_prev, ph_max]
+                     ! ph_new = (ph_prev + ph_max)/2d0
+                     ! In terms of [H]_new:
+                     ! [H]_new = 10**(-ph_new)
+                     !         = 10**(-(ph_prev + ph_max)/2d0)
+                     !         = SQRT(10**(-(ph_prev + phmax)))
+                     !         = SQRT([H]_old*10**(-ph_max))
+                     !         = SQRT([H]_old * zh_min)
+                     zh                = SQRT(zh_prev * zh_min(ji,jj,jk))
+                     zh_lnfactor       = (zh - zh_prev)/zh_prev ! Required to test convergence below
+                  ENDIF
+
+                  IF( zh > zh_max(ji,jj,jk) ) THEN
+                     ! if [H]_new > [H]_max
+                     ! i.e., if ph_new < ph_min, then
+                     ! take one bisection step on [ph_min, ph_prev]
+                     ! ph_new = (ph_prev + ph_min)/2d0
+                     ! In terms of [H]_new:
+                     ! [H]_new = 10**(-ph_new)
+                     !         = 10**(-(ph_prev + ph_min)/2d0)
+                     !         = SQRT(10**(-(ph_prev + ph_min)))
+                     !         = SQRT([H]_old*10**(-ph_min))
+                     !         = SQRT([H]_old * zhmax)
+                     zh                = SQRT(zh_prev * zh_max(ji,jj,jk))
+                     zh_lnfactor       = (zh - zh_prev)/zh_prev ! Required to test convergence below
+                  ENDIF
+               ENDIF
+
+               zeqn_absmin(ji,jj,jk) = MIN( ABS(zeqn), zeqn_absmin(ji,jj,jk))
+
+               ! Stop iterations once |\delta{[H]}/[H]| < rdel
+               ! <=> |(zh - zh_prev)/zh_prev| = |EXP(-zeqn/(zdeqndh*zh_prev)) -1| < rdel
+               ! |EXP(-zeqn/(zdeqndh*zh_prev)) -1| ~ |zeqn/(zdeqndh*zh_prev)|
+
+               ! Alternatively:
+               ! |\Delta pH| = |zeqn/(zdeqndh*zh_prev*LOG(10))|
+               !             ~ 1/LOG(10) * |\Delta [H]|/[H]
+               !             < 1/LOG(10) * rdel
+
+               ! Hence |zeqn/(zdeqndh*zh)| < rdel
+
+               ! rdel <-- pp_rdel_ah_target
+               l_exitnow = (ABS(zh_lnfactor) < pp_rdel_ah_target)
+
+               IF(l_exitnow) THEN 
+                  rmask(ji,jj,jk) = 0.
+               ENDIF
+
+               zhi(ji,jj,jk) =  zh
+
+               IF(jn >= jp_maxniter_atgen) THEN
+                  zhi(ji,jj,jk) = -1._wp
+               ENDIF
+
+            ENDIF
+         END DO
+      END DO
+   END DO
+   END DO
+   !
+   CALL wrk_dealloc( jpi, jpj, jpk, zalknw_inf, zalknw_sup, rmask )
+   CALL wrk_dealloc( jpi, jpj, jpk, zh_min, zh_max, zeqn_absmin )
+
+
+   IF( nn_timing == 1 )  CALL timing_stop('p4z_che_solve_hi')
+
+
+   END SUBROUTINE p4z_che_solve_hi
 
    INTEGER FUNCTION p4z_che_alloc()
@@ -336,6 +817,23 @@
       !!                     ***  ROUTINE p4z_che_alloc  ***
       !!----------------------------------------------------------------------
-      ALLOCATE( sio3eq(jpi,jpj,jpk), fekeq(jpi,jpj,jpk), chemc(jpi,jpj,3), chemo2(jpi,jpj,jpk),   &
-      &         tempis(jpi,jpj,jpk), STAT=p4z_che_alloc )
+      INTEGER ::   ierr(3)        ! Local variables
+      !!----------------------------------------------------------------------
+
+      ierr(:) = 0
+
+      ALLOCATE( sio3eq(jpi,jpj,jpk), fekeq(jpi,jpj,jpk), chemc(jpi,jpj,3), chemo2(jpi,jpj,jpk), STAT=ierr(1) )
+
+      ALLOCATE( akb3(jpi,jpj,jpk)     , tempis(jpi, jpj, jpk),       &
+         &      akw3(jpi,jpj,jpk)     , borat (jpi,jpj,jpk)  ,       &
+         &      aks3(jpi,jpj,jpk)     , akf3(jpi,jpj,jpk)    ,       &
+         &      ak1p3(jpi,jpj,jpk)    , ak2p3(jpi,jpj,jpk)   ,       &
+         &      ak3p3(jpi,jpj,jpk)    , aksi3(jpi,jpj,jpk)   ,       &
+         &      fluorid(jpi,jpj,jpk)  , sulfat(jpi,jpj,jpk)  ,       &
+         &      salinprac(jpi,jpj,jpk),                 STAT=ierr(2) )
+
+      ALLOCATE( fesol(jpi,jpj,jpk,5), STAT=ierr(3) )
+
+      !* Variable for chemistry of the CO2 cycle
+      p4z_che_alloc = MAXVAL( ierr )
       !
       IF( p4z_che_alloc /= 0 )   CALL ctl_warn('p4z_che_alloc : failed to allocate arrays.')
@@ -355,3 +853,3 @@
 
    !!======================================================================
-END MODULE p4zche
+END MODULE  p4zche
Index: branches/2015/nemo_v3_6_STABLE/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zflx.F90
===================================================================
--- branches/2015/nemo_v3_6_STABLE/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zflx.F90	(revision 7522)
+++ branches/2015/nemo_v3_6_STABLE/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zflx.F90	(revision 7607)
@@ -87,5 +87,5 @@
       REAL(wp) ::   zfld, zflu, zfld16, zflu16, zfact
       REAL(wp) ::   zvapsw, zsal, zfco2, zxc2, xCO2approx, ztkel, zfugcoeff
-      REAL(wp) ::   zph, zah2, zbot, zdic, zalk, zsch_o2, zalka, zsch_co2
+      REAL(wp) ::   zph, zdic, zsch_o2, zsch_co2
       REAL(wp) ::   zyr_dec, zdco2dt
       CHARACTER (len=25) :: charout
@@ -122,31 +122,14 @@
 #endif
 
-      DO jm = 1, 10
-!CDIR NOVERRCHK
-         DO jj = 1, jpj
-!CDIR NOVERRCHK
-            DO ji = 1, jpi
-
-               ! DUMMY VARIABLES FOR DIC, H+, AND BORATE
-               zbot  = borat(ji,jj,1)
-               zfact = rhop(ji,jj,1) / 1000. + rtrn
-               zdic  = trb(ji,jj,1,jpdic) / zfact
-               zph   = MAX( hi(ji,jj,1), 1.e-10 ) / zfact
-               zalka = trb(ji,jj,1,jptal) / zfact
-
-               ! CALCULATE [ALK]([CO3--], [HCO3-])
-               zalk  = zalka - (  akw3(ji,jj,1) / zph - zph / aphscale(ji,jj,1)    &
-               &       + zbot / ( 1.+ zph / akb3(ji,jj,1) )  )
-
-               ! CALCULATE [H+] AND [H2CO3]
-               zah2   = SQRT(  (zdic-zalk)**2 + 4.* ( zalk * ak23(ji,jj,1)   &
-                  &                                        / ak13(ji,jj,1) ) * ( 2.* zdic - zalk )  )
-               zah2   = 0.5 * ak13(ji,jj,1) / zalk * ( ( zdic - zalk ) + zah2 )
-               zh2co3(ji,jj) = ( 2.* zdic - zalk ) / ( 2.+ ak13(ji,jj,1) / zah2 ) * zfact
-               hi(ji,jj,1)   = zah2 * zfact
-            END DO
+      DO jj = 1, jpj
+         DO ji = 1, jpi
+            ! DUMMY VARIABLES FOR DIC, H+, AND BORATE
+            zfact = rhop(ji,jj,1) / 1000. + rtrn
+            zdic  = trb(ji,jj,1,jpdic)
+            zph   = MAX( hi(ji,jj,1), 1.e-10 ) / zfact
+            ! CALCULATE [H2CO3]
+            zh2co3(ji,jj) = zdic/(1. + ak13(ji,jj,1)/zph + ak13(ji,jj,1)*ak23(ji,jj,1)/zph**2)
          END DO
       END DO
-
 
       ! --------------
@@ -254,4 +237,9 @@
       ENDIF
       !
+#if defined key_cpl_carbon_cycle
+      ! change units for carbon cycle coupling
+      oce_co2(:,:) = oce_co2(:,:) / e1e2t(:,:) * rfact2r ! in molC/m2/s
+#endif
+      !
       CALL wrk_dealloc( jpi, jpj, zkgco2, zkgo2, zh2co3, zoflx, zpco2atm )
       !
Index: branches/2015/nemo_v3_6_STABLE/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zlys.F90
===================================================================
--- branches/2015/nemo_v3_6_STABLE/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zlys.F90	(revision 7522)
+++ branches/2015/nemo_v3_6_STABLE/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zlys.F90	(revision 7607)
@@ -23,4 +23,5 @@
    USE sms_pisces      !  PISCES Source Minus Sink variables
    USE prtctl_trc      !  print control for debugging
+   USE p4zche          !  Chemical model
    USE iom             !  I/O manager
 
@@ -60,49 +61,33 @@
       !
       INTEGER, INTENT(in) ::   kt, knt ! ocean time step
-      INTEGER  ::   ji, jj, jk, jn
-      REAL(wp) ::   zalk, zdic, zph, zah2
-      REAL(wp) ::   zdispot, zfact, zcalcon, zalka, zaldi
+      INTEGER  ::   ji, jj, jk
+      REAL(wp) ::   zdispot, zfact, zcalcon
       REAL(wp) ::   zomegaca, zexcess, zexcess0
       CHARACTER (len=25) :: charout
-      REAL(wp), POINTER, DIMENSION(:,:,:) :: zco3, zco3sat, zcaldiss   
+      REAL(wp), POINTER, DIMENSION(:,:,:) :: zco3, zco3sat, zcaldiss, zhinit, zhi   
       !!---------------------------------------------------------------------
       !
       IF( nn_timing == 1 )  CALL timing_start('p4z_lys')
       !
-      CALL wrk_alloc( jpi, jpj, jpk, zco3, zco3sat, zcaldiss )
+      CALL wrk_alloc( jpi, jpj, jpk, zco3, zco3sat, zcaldiss, zhinit, zhi )
       !
       zco3    (:,:,:) = 0.
       zcaldiss(:,:,:) = 0.
+      zhinit(:,:,:)   = hi(:,:,:) * 1000. / ( rhop(:,:,:) + rtrn )
       !     -------------------------------------------
       !     COMPUTE [CO3--] and [H+] CONCENTRATIONS
       !     -------------------------------------------
-      
-      DO jn = 1, 5                               !  BEGIN OF ITERATION
-         !
-!CDIR NOVERRCHK
-         DO jk = 1, jpkm1
-!CDIR NOVERRCHK
-            DO jj = 1, jpj
-!CDIR NOVERRCHK
-               DO ji = 1, jpi
-                  zfact = rhop(ji,jj,jk) / 1000. + rtrn
-                  zph  = hi(ji,jj,jk) * tmask(ji,jj,jk) / zfact + ( 1.-tmask(ji,jj,jk) ) * 1.e-9 ! [H+]
-                  zdic  = trb(ji,jj,jk,jpdic) / zfact
-                  zalka = trb(ji,jj,jk,jptal) / zfact
-                  ! CALCULATE [ALK]([CO3--], [HCO3-])
-                  zalk  = zalka - ( akw3(ji,jj,jk) / zph - zph / ( aphscale(ji,jj,jk) + rtrn )  &
-                  &       + borat(ji,jj,jk) / ( 1. + zph / akb3(ji,jj,jk) ) )
-                  ! CALCULATE [H+] and [CO3--]
-                  zaldi = zdic - zalk
-                  zah2  = SQRT( zaldi * zaldi + 4.* ( zalk * ak23(ji,jj,jk) / ak13(ji,jj,jk) ) * ( zdic + zaldi ) )
-                  zah2  = 0.5 * ak13(ji,jj,jk) / zalk * ( zaldi + zah2 )
-                  !
-                  zco3(ji,jj,jk) = zalk / ( 2. + zah2 / ak23(ji,jj,jk) ) * zfact
-                  hi(ji,jj,jk)   = zah2 * zfact
-               END DO
+
+      CALL p4z_che_solve_hi( zhinit, zhi )
+
+      DO jk = 1, jpkm1
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               zco3(ji,jj,jk) = trb(ji,jj,jk,jpdic) * ak13(ji,jj,jk) * ak23(ji,jj,jk) / (zhi(ji,jj,jk)**2   &
+               &                + ak13(ji,jj,jk) * zhi(ji,jj,jk) + ak13(ji,jj,jk) * ak23(ji,jj,jk) + rtrn )
+               hi(ji,jj,jk)   = zhi(ji,jj,jk) * rhop(ji,jj,jk) / 1000.
             END DO
          END DO
-         !
-      END DO 
+      END DO
 
       !     ---------------------------------------------------------
@@ -138,5 +123,4 @@
               !       AND [SUM(CO2)] DUE TO CACO3 DISSOLUTION/PRECIPITATION
               zcaldiss(ji,jj,jk)  = zdispot * rfact2 / rmtss ! calcite dissolution
-              zco3(ji,jj,jk)      = zco3(ji,jj,jk) + zcaldiss(ji,jj,jk)
               !
               tra(ji,jj,jk,jptal) = tra(ji,jj,jk,jptal) + 2. * zcaldiss(ji,jj,jk)
@@ -167,5 +151,5 @@
       ENDIF
       !
-      CALL wrk_dealloc( jpi, jpj, jpk, zco3, zco3sat, zcaldiss )
+      CALL wrk_dealloc( jpi, jpj, jpk, zco3, zcaldiss, zhinit, zhi, zco3sat )
       !
       IF( nn_timing == 1 )  CALL timing_stop('p4z_lys')
@@ -186,9 +170,5 @@
       !!
       !!----------------------------------------------------------------------
-      INTEGER  ::  ji, jj, jk
       INTEGER  ::  ios                 ! Local integer output status for namelist read
-      REAL(wp) ::  zcaralk, zbicarb, zco3
-      REAL(wp) ::  ztmas, ztmas1
-
       NAMELIST/nampiscal/ kdca, nca
       !!----------------------------------------------------------------------
Index: branches/2015/nemo_v3_6_STABLE/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zsms.F90
===================================================================
--- branches/2015/nemo_v3_6_STABLE/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zsms.F90	(revision 7522)
+++ branches/2015/nemo_v3_6_STABLE/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zsms.F90	(revision 7607)
@@ -85,5 +85,5 @@
         CALL p4z_che                              ! initialize the chemical constants
         !
-        IF( .NOT. ln_rsttr ) THEN  ;   CALL p4z_ph_ini   !  set PH at kt=nit000 
+        IF( .NOT. ln_rsttr ) THEN  ;   CALL p4z_che_ahini( hi )   !  set PH at kt=nit000
         ELSE                       ;   CALL p4z_rst( nittrc000, 'READ' )  !* read or initialize all required fields 
         ENDIF
@@ -310,32 +310,4 @@
    END SUBROUTINE p4z_sms_init
 
-   SUBROUTINE p4z_ph_ini
-      !!---------------------------------------------------------------------
-      !!                   ***  ROUTINE p4z_ini_ph  ***
-      !!
-      !!  ** Purpose : Initialization of chemical variables of the carbon cycle
-      !!---------------------------------------------------------------------
-      INTEGER  ::  ji, jj, jk
-      REAL(wp) ::  zcaralk, zbicarb, zco3
-      REAL(wp) ::  ztmas, ztmas1
-      !!---------------------------------------------------------------------
-
-      ! Set PH from  total alkalinity, borat (???), akb3 (???) and ak23 (???)
-      ! --------------------------------------------------------
-      DO jk = 1, jpk
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-               ztmas   = tmask(ji,jj,jk)
-               ztmas1  = 1. - tmask(ji,jj,jk)
-               zcaralk = trb(ji,jj,jk,jptal) - borat(ji,jj,jk) / (  1. + 1.E-8 / ( rtrn + akb3(ji,jj,jk) )  )
-               zco3    = ( zcaralk - trb(ji,jj,jk,jpdic) ) * ztmas + 0.5e-3 * ztmas1
-               zbicarb = ( 2. * trb(ji,jj,jk,jpdic) - zcaralk )
-               hi(ji,jj,jk) = ( ak23(ji,jj,jk) * zbicarb / zco3 ) * ztmas + 1.e-9 * ztmas1
-            END DO
-         END DO
-     END DO
-     !
-   END SUBROUTINE p4z_ph_ini
-
    SUBROUTINE p4z_rst( kt, cdrw )
       !!---------------------------------------------------------------------
@@ -350,8 +322,4 @@
       INTEGER         , INTENT(in) ::   kt         ! ocean time-step
       CHARACTER(len=*), INTENT(in) ::   cdrw       ! "READ"/"WRITE" flag
-      !
-      INTEGER  ::  ji, jj, jk
-      REAL(wp) ::  zcaralk, zbicarb, zco3
-      REAL(wp) ::  ztmas, ztmas1
       !!---------------------------------------------------------------------
 
@@ -365,6 +333,5 @@
             CALL iom_get( numrtr, jpdom_autoglo, 'PH' , hi(:,:,:)  )
          ELSE
-!            hi(:,:,:) = 1.e-9 
-            CALL p4z_ph_ini
+            CALL p4z_che_ahini( hi )
          ENDIF
          CALL iom_get( numrtr, jpdom_autoglo, 'Silicalim', xksi(:,:) )
Index: branches/2015/nemo_v3_6_STABLE/NEMOGCM/NEMO/TOP_SRC/PISCES/sms_pisces.F90
===================================================================
--- branches/2015/nemo_v3_6_STABLE/NEMOGCM/NEMO/TOP_SRC/PISCES/sms_pisces.F90	(revision 7522)
+++ branches/2015/nemo_v3_6_STABLE/NEMOGCM/NEMO/TOP_SRC/PISCES/sms_pisces.F90	(revision 7607)
@@ -93,10 +93,7 @@
 
    !!* Variable for chemistry of the CO2 cycle
-   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   akb3       !: ???
    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   ak13       !: ???
    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   ak23       !: ???
    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   aksp       !: ???
-   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   akw3       !: ???
-   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   borat      !: ???
    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   hi         !: ???
    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   excess     !: ???
@@ -153,7 +150,6 @@
 
       !* Variable for chemistry of the CO2 cycle
-      ALLOCATE( akb3(jpi,jpj,jpk)    , ak13  (jpi,jpj,jpk) ,       &
+      ALLOCATE( ak13  (jpi,jpj,jpk) ,                              &
          &      ak23(jpi,jpj,jpk)    , aksp  (jpi,jpj,jpk) ,       &
-         &      akw3(jpi,jpj,jpk)    , borat (jpi,jpj,jpk) ,       &
          &      hi  (jpi,jpj,jpk)    , excess(jpi,jpj,jpk) ,       &
          &      aphscale(jpi,jpj,jpk),                           STAT=ierr(4) )
Index: branches/2015/nemo_v3_6_STABLE/NEMOGCM/NEMO/TOP_SRC/trcstp.F90
===================================================================
--- branches/2015/nemo_v3_6_STABLE/NEMOGCM/NEMO/TOP_SRC/trcstp.F90	(revision 7522)
+++ branches/2015/nemo_v3_6_STABLE/NEMOGCM/NEMO/TOP_SRC/trcstp.F90	(revision 7607)
@@ -131,5 +131,5 @@
       INTEGER, INTENT(in) ::   kt
       INTEGER  :: jn
-      REAL(wp) :: zkt
+      REAL(wp) :: zkt, zrec
       CHARACTER(len=1)               ::   cl1                      ! 1 character
       CHARACTER(len=2)               ::   cl2                      ! 2 characters
@@ -153,20 +153,29 @@
          !
          !                                            !* Restart: read in restart file
-         IF( ln_rsttr .AND. iom_varid( numrtr, 'qsr_mean' , ldstop = .FALSE. ) > 0 .AND. &
-                            iom_varid( numrtr, 'qsr_arr_1', ldstop = .FALSE. ) > 0 .AND. &
-                            iom_varid( numrtr, 'ktdcy'    , ldstop = .FALSE. ) > 0 ) THEN 
-            CALL iom_get( numrtr, 'ktdcy', zkt )   !  A mean of qsr
+         IF(  ln_rsttr .AND. iom_varid( numrtr, 'qsr_mean' , ldstop = .FALSE. ) > 0  &
+           &           .AND. iom_varid( numrtr, 'qsr_arr_1', ldstop = .FALSE. ) > 0  &
+           &           .AND. iom_varid( numrtr, 'ktdcy'    , ldstop = .FALSE. ) > 0  &
+           &           .AND. iom_varid( numrtr, 'nrdcy'    , ldstop = .FALSE. ) > 0  ) THEN
+
+            CALL iom_get( numrtr, 'ktdcy', zkt )  
             rsecfst = INT( zkt ) * rdttrc(1)
             IF(lwp) WRITE(numout,*) 'trc_qsr_mean:   qsr_mean read in the restart file at time-step rsecfst =', rsecfst, ' s '
             CALL iom_get( numrtr, jpdom_autoglo, 'qsr_mean', qsr_mean )   !  A mean of qsr
-            DO jn = 1, nb_rec_per_day 
-             IF( jn <= 9 )  THEN
-               WRITE(cl1,'(i1)') jn
-               CALL iom_get( numrtr, jpdom_autoglo, 'qsr_arr_'//cl1, qsr_arr(:,:,jn) )   !  A mean of qsr
-             ELSE
-               WRITE(cl2,'(i2.2)') jn
-               CALL iom_get( numrtr, jpdom_autoglo, 'qsr_arr_'//cl2, qsr_arr(:,:,jn) )   !  A mean of qsr
-             ENDIF
-           ENDDO
+            CALL iom_get( numrtr, 'nrdcy', zrec )   !  Number of record per days
+            IF( INT( zrec ) == nb_rec_per_day ) THEN
+               DO jn = 1, nb_rec_per_day 
+                  IF( jn <= 9 )  THEN
+                    WRITE(cl1,'(i1)') jn
+                    CALL iom_get( numrtr, jpdom_autoglo, 'qsr_arr_'//cl1, qsr_arr(:,:,jn) )   !  A mean of qsr
+                  ELSE
+                    WRITE(cl2,'(i2.2)') jn
+                    CALL iom_get( numrtr, jpdom_autoglo, 'qsr_arr_'//cl2, qsr_arr(:,:,jn) )   !  A mean of qsr
+                  ENDIF
+              ENDDO
+            ELSE
+               DO jn = 1, nb_rec_per_day
+                  qsr_arr(:,:,jn) = qsr_mean(:,:)
+               ENDDO
+            ENDIF
          ELSE                                         !* no restart: set from nit000 values
             IF(lwp) WRITE(numout,*) 'trc_qsr_mean:   qsr_mean set to nit000 values'
@@ -185,5 +194,5 @@
       llnew   = ( rseclast - rsecfst ) .ge.  rdt_sampl    !   new shortwave to store
       IF( llnew ) THEN
-          IF( lwp ) WRITE(numout,*) ' New shortwave to sample for TOP at time kt = ', kt, &
+          IF( lwp .AND. kt < nittrc000 + 100 ) WRITE(numout,*) ' New shortwave to sample for TOP at time kt = ', kt, &
              &                      ' time = ', rseclast/3600.,'hours '
           rsecfst = rseclast
@@ -199,6 +208,8 @@
          IF(lwp) WRITE(numout,*) 'trc_mean_qsr : write qsr_mean in restart file  kt =', kt
          IF(lwp) WRITE(numout,*) '~~~~~~~'
-         zkt = REAL( kt, wp )
-         CALL iom_rstput( kt, nitrst, numrtw, 'ktdcy', zkt )
+         zkt  = REAL( kt, wp )
+         zrec = REAL( nb_rec_per_day, wp )
+         CALL iom_rstput( kt, nitrst, numrtw, 'ktdcy', zkt  )
+         CALL iom_rstput( kt, nitrst, numrtw, 'nrdcy', zrec )
           DO jn = 1, nb_rec_per_day 
              IF( jn <= 9 )  THEN
