Index: branches/2015/nemo_v3_6_STABLE/NEMOGCM/CONFIG/ORCA2_LIM3/EXP00/namelist_ice_cfg
===================================================================
--- branches/2015/nemo_v3_6_STABLE/NEMOGCM/CONFIG/ORCA2_LIM3/EXP00/namelist_ice_cfg	(revision 7564)
+++ branches/2015/nemo_v3_6_STABLE/NEMOGCM/CONFIG/ORCA2_LIM3/EXP00/namelist_ice_cfg	(revision 7607)
@@ -1,4 +1,4 @@
 !!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
-!! NEMO/LIM-2 : Ice configuration namelist. Overwrites SHARED/namelist_ice_lim2_ref
+!! NEMO/LIM-3 : Ice configuration namelist. Overwrites SHARED/namelist_ice_lim3_ref
 !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
 
@@ -15,4 +15,9 @@
 !-----------------------------------------------------------------------
 /
+!------------------------------------------------------------------------------
+&namicehdf     !   Ice horizontal diffusion
+!------------------------------------------------------------------------------
+   nn_ahi0        =  2          !  horizontal diffusivity computation
+/
 !-----------------------------------------------------------------------
 &namicethd     !   ice thermodynamic
@@ -24,5 +29,5 @@
 /
 !-----------------------------------------------------------------------
-&namiceitdme   !   parameters for mechanical redistribution of ice 
+&namiceitdme   !   parameters for mechanical redistribution of ice
 !-----------------------------------------------------------------------
 /
@@ -31,2 +36,6 @@
 !-----------------------------------------------------------------------
 /
+!------------------------------------------------------------------------------
+&namiceitd     !   Ice discretization
+!------------------------------------------------------------------------------
+/
Index: branches/2015/nemo_v3_6_STABLE/NEMOGCM/CONFIG/SHARED/field_def.xml
===================================================================
--- branches/2015/nemo_v3_6_STABLE/NEMOGCM/CONFIG/SHARED/field_def.xml	(revision 7564)
+++ branches/2015/nemo_v3_6_STABLE/NEMOGCM/CONFIG/SHARED/field_def.xml	(revision 7607)
@@ -225,4 +225,6 @@
          <field id="icealb_cea"   long_name="Ice albedo (cell average)"                 standard_name="sea_ice_albedo"                            unit="1"        />
          <field id="calving_cea"  long_name="Calving"                                   standard_name="water_flux_into_sea_water_from_icebergs"   unit="kg/m2/s"  />
+         <field id="iceberg_cea"  long_name="Iceberg"                                   standard_name="water_flux_into_sea_water_from_icebergs"   unit="kg/m2/s"  />
+         <field id="iceshelf_cea"  long_name="Iceshelf"                                   standard_name="water_flux_into_sea_water_from_iceshelf"   unit="kg/m2/s"  />
 
          <!-- available if key_oasis3 + conservative method -->
Index: branches/2015/nemo_v3_6_STABLE/NEMOGCM/CONFIG/SHARED/namelist_ice_lim3_ref
===================================================================
--- branches/2015/nemo_v3_6_STABLE/NEMOGCM/CONFIG/SHARED/namelist_ice_lim3_ref	(revision 7564)
+++ branches/2015/nemo_v3_6_STABLE/NEMOGCM/CONFIG/SHARED/namelist_ice_lim3_ref	(revision 7607)
@@ -69,14 +69,15 @@
    rn_relast      =    0.333       !  ratio of elastic timescale to ice time step: Telast = dt_ice * rn_relast 
                                    !     advised value: 1/3 (rn_nevp=120) or 1/9 (rn_nevp=300)
-   nn_ahi0        =    2           !  horizontal diffusivity computation
-                                   !     0: use rn_ahi0_ref
-                                   !     1: use rn_ahi0_ref x mean grid cell length / ( 2deg mean grid cell length )
-                                   !     2: use rn_ahi0_ref x grid cell length      / ( 2deg mean grid cell length )
-   rn_ahi0_ref    = 350.0          !  horizontal sea ice diffusivity (m2/s) 
-                                   !     if nn_ahi0 > 0, rn_ahi0_ref is the reference value at a nominal 2 deg resolution
 /
 !------------------------------------------------------------------------------
 &namicehdf     !   Ice horizontal diffusion
 !------------------------------------------------------------------------------
+   nn_ahi0        =    -1          !  horizontal diffusivity computation
+                                   !    -1: no diffusion (bypass limhdf)
+                                   !     0: use rn_ahi0_ref
+                                   !     1: use rn_ahi0_ref x mean grid cell length / ( 2deg mean grid cell length )
+                                   !     2: use rn_ahi0_ref x grid cell length      / ( 2deg mean grid cell length )
+   rn_ahi0_ref    = 350.0          !  horizontal sea ice diffusivity (m2/s)
+                                   !     if nn_ahi0 > 0, rn_ahi0_ref is the reference value at a nominal 2 deg resolution
    nn_convfrq     = 5              !  convergence check frequency of the Crant-Nicholson scheme (perf. optimization)
 /
@@ -98,4 +99,6 @@
                                    !     0: k = k0 + beta.S/T (Untersteiner, 1964)
                                    !     1: k = k0 + beta1.S/T - beta2.T (Pringle et al., 2007)
+   rn_cdsn     = 0.31              !  thermal conductivity of the snow (0.31 W/m/K, Maykut and Untersteiner, 1971)
+                                   !  Obs: 0.1-0.5 (Lecomte et al, JAMES 2013)
    nn_monocat  = 0                 !  virtual ITD mono-category parameterizations (1, jpl = 1 only) or not (0)
                                    !     2: simple piling instead of ridging --- temporary option
Index: branches/2015/nemo_v3_6_STABLE/NEMOGCM/CONFIG/SHARED/namelist_ref
===================================================================
--- branches/2015/nemo_v3_6_STABLE/NEMOGCM/CONFIG/SHARED/namelist_ref	(revision 7564)
+++ branches/2015/nemo_v3_6_STABLE/NEMOGCM/CONFIG/SHARED/namelist_ref	(revision 7607)
@@ -500,11 +500,12 @@
 &namsbc_alb    !   albedo parameters
 !-----------------------------------------------------------------------
-   nn_ice_alb  =    0   !  parameterization of ice/snow albedo
-                        !     0: Shine & Henderson-Sellers (JGR 1985)
-                        !     1: "home made" based on Brandt et al. (J. Climate 2005)
-                        !                         and Grenfell & Perovich (JGR 2004)
-   rn_albice   =  0.53  !  albedo of bare puddled ice (values from 0.49 to 0.58)
-                        !     0.53 (default) => if nn_ice_alb=0
-                        !     0.50 (default) => if nn_ice_alb=1
+   nn_ice_alb   =    1   !  parameterization of ice/snow albedo
+                         !     0: Shine & Henderson-Sellers (JGR 1985), giving clear-sky albedo
+                         !     1: "home made" based on Brandt et al. (JClim 2005) and Grenfell & Perovich (JGR 2004),
+                         !        giving cloud-sky albedo
+   rn_alb_sdry  =  0.85  !  dry snow albedo         : 0.80 (nn_ice_alb = 0); 0.85 (nn_ice_alb = 1); obs 0.85-0.87 (cloud-sky)
+   rn_alb_smlt  =  0.75  !  melting snow albedo     : 0.65 ( '' )          ; 0.75 ( '' )          ; obs 0.72-0.82 ( '' )
+   rn_alb_idry  =  0.60  !  dry ice albedo          : 0.72 ( '' )          ; 0.60 ( '' )          ; obs 0.54-0.65 ( '' )
+   rn_alb_imlt  =  0.50  !  bare puddled ice albedo : 0.53 ( '' )          ; 0.50 ( '' )          ; obs 0.49-0.58 ( '' )
 /
 !-----------------------------------------------------------------------
Index: branches/2015/nemo_v3_6_STABLE/NEMOGCM/NEMO/LIM_SRC_3/ice.F90
===================================================================
--- branches/2015/nemo_v3_6_STABLE/NEMOGCM/NEMO/LIM_SRC_3/ice.F90	(revision 7564)
+++ branches/2015/nemo_v3_6_STABLE/NEMOGCM/NEMO/LIM_SRC_3/ice.F90	(revision 7607)
@@ -212,4 +212,5 @@
    REAL(wp), PUBLIC ::   rn_betas         !: coef. for partitioning of snowfall between leads and sea ice
    REAL(wp), PUBLIC ::   rn_kappa_i       !: coef. for the extinction of radiation Grenfell et al. (2006) [1/m]
+   REAL(wp), PUBLIC ::   rn_cdsn          !: thermal conductivity of the snow [W/m/K]
    REAL(wp), PUBLIC ::   nn_conv_dif      !: maximal number of iterations for heat diffusion
    REAL(wp), PUBLIC ::   rn_terr_dif      !: maximal tolerated error (C) for heat diffusion
@@ -320,4 +321,5 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   o_i     !: Sea-Ice Age (days)
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   oa_i    !: Sea-Ice Age times ice area (days)
+
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   bv_i    !: brine volume
 
@@ -406,4 +408,5 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   diag_vice     !: ice volume variation   [m/s] 
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   diag_vsnw     !: snw volume variation   [m/s] 
+
    !
    !!----------------------------------------------------------------------
@@ -463,5 +466,5 @@
          &      et_i (jpi,jpj) , et_s (jpi,jpj) , tm_i (jpi,jpj) , bvm_i(jpi,jpj) ,     &
          &      smt_i(jpi,jpj) , tm_su(jpi,jpj) , htm_i(jpi,jpj) , htm_s(jpi,jpj) ,     &
-         &      om_i (jpi,jpj)                              , STAT=ierr(ii) )
+         &      om_i (jpi,jpj) , STAT=ierr(ii) )
       ii = ii + 1
       ALLOCATE( t_s(jpi,jpj,nlay_s,jpl) , e_s(jpi,jpj,nlay_s,jpl) , STAT=ierr(ii) )
Index: branches/2015/nemo_v3_6_STABLE/NEMOGCM/NEMO/LIM_SRC_3/limdyn.F90
===================================================================
--- branches/2015/nemo_v3_6_STABLE/NEMOGCM/NEMO/LIM_SRC_3/limdyn.F90	(revision 7564)
+++ branches/2015/nemo_v3_6_STABLE/NEMOGCM/NEMO/LIM_SRC_3/limdyn.F90	(revision 7607)
@@ -244,7 +244,5 @@
       INTEGER  ::   ios                 ! Local integer output status for namelist read
       NAMELIST/namicedyn/ nn_icestr, ln_icestr_bvf, rn_pe_rdg, rn_pstar, rn_crhg, rn_cio,  rn_creepl, rn_ecc, &
-         &                nn_nevp, rn_relast, nn_ahi0, rn_ahi0_ref
-      INTEGER  ::   ji, jj
-      REAL(wp) ::   za00, zd_max
+         &                nn_nevp, rn_relast
       !!-------------------------------------------------------------------
 
@@ -272,6 +270,4 @@
          WRITE(numout,*) '   number of iterations for subcycling                  nn_nevp       = ', nn_nevp
          WRITE(numout,*) '   ratio of elastic timescale over ice time step        rn_relast     = ', rn_relast
-         WRITE(numout,*) '   horizontal diffusivity calculation                   nn_ahi0       = ', nn_ahi0
-         WRITE(numout,*) '   horizontal diffusivity coeff. (orca2 grid)           rn_ahi0_ref   = ', rn_ahi0_ref
       ENDIF
       !
@@ -279,47 +275,4 @@
       rhoco  = rau0  * rn_cio
       !
-      !  Diffusion coefficients
-      SELECT CASE( nn_ahi0 )
-
-      CASE( 0 )
-         ahiu(:,:) = rn_ahi0_ref
-         ahiv(:,:) = rn_ahi0_ref
-
-         IF(lwp) WRITE(numout,*) ''
-         IF(lwp) WRITE(numout,*) '   laplacian operator: ahim constant = rn_ahi0_ref'
-
-      CASE( 1 ) 
-
-         zd_max = MAX( MAXVAL( e1t(:,:) ), MAXVAL( e2t(:,:) ) )
-         IF( lk_mpp )   CALL mpp_max( zd_max )          ! max over the global domain
-         
-         ahiu(:,:) = rn_ahi0_ref * zd_max * 1.e-05_wp   ! 1.e05 = 100km = max grid space at 60° latitude in orca2
-                                                        !                    (60° = min latitude for ice cover)  
-         ahiv(:,:) = rn_ahi0_ref * zd_max * 1.e-05_wp
-
-         IF(lwp) WRITE(numout,*) ''
-         IF(lwp) WRITE(numout,*) '   laplacian operator: ahim proportional to max of e1 e2 over the domain (', zd_max, ')'
-         IF(lwp) WRITE(numout,*) '   value for ahim = ', rn_ahi0_ref * zd_max * 1.e-05_wp 
-         
-      CASE( 2 ) 
-
-         zd_max = MAX( MAXVAL( e1t(:,:) ), MAXVAL( e2t(:,:) ) )
-         IF( lk_mpp )   CALL mpp_max( zd_max )   ! max over the global domain
-         
-         za00 = rn_ahi0_ref * 1.e-05_wp          ! 1.e05 = 100km = max grid space at 60° latitude in orca2
-                                                 !                    (60° = min latitude for ice cover)  
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-               ahiu(ji,jj) = za00 * MAX( e1t(ji,jj), e2t(ji,jj) ) * umask(ji,jj,1)
-               ahiv(ji,jj) = za00 * MAX( e1f(ji,jj), e2f(ji,jj) ) * vmask(ji,jj,1)
-            END DO
-         END DO
-         !
-         IF(lwp) WRITE(numout,*) ''
-         IF(lwp) WRITE(numout,*) '   laplacian operator: ahim proportional to e1'
-         IF(lwp) WRITE(numout,*) '   maximum grid-spacing = ', zd_max, ' maximum value for ahim = ', za00*zd_max
-         
-      END SELECT
-
    END SUBROUTINE lim_dyn_init
 
Index: branches/2015/nemo_v3_6_STABLE/NEMOGCM/NEMO/LIM_SRC_3/limhdf.F90
===================================================================
--- branches/2015/nemo_v3_6_STABLE/NEMOGCM/NEMO/LIM_SRC_3/limhdf.F90	(revision 7564)
+++ branches/2015/nemo_v3_6_STABLE/NEMOGCM/NEMO/LIM_SRC_3/limhdf.F90	(revision 7607)
@@ -233,12 +233,8 @@
       !!-------------------------------------------------------------------
       INTEGER  ::   ios                 ! Local integer output status for namelist read
-      NAMELIST/namicehdf/ nn_convfrq 
-      !!-------------------------------------------------------------------
-      !
-      IF(lwp) THEN
-         WRITE(numout,*)
-         WRITE(numout,*) 'lim_hdf : Ice horizontal diffusion'
-         WRITE(numout,*) '~~~~~~~'
-      ENDIF
+      NAMELIST/namicehdf/  nn_ahi0, rn_ahi0_ref, nn_convfrq 
+      INTEGER  ::   ji, jj
+      REAL(wp) ::   za00, zd_max
+      !!-------------------------------------------------------------------
       !
       REWIND( numnam_ice_ref )              ! Namelist namicehdf in reference namelist : Ice horizontal diffusion
@@ -253,7 +249,60 @@
       IF(lwp) THEN                          ! control print
          WRITE(numout,*)
-         WRITE(numout,*)'   Namelist of ice parameters for ice horizontal diffusion computation '
-         WRITE(numout,*)'      convergence check frequency of the Crant-Nicholson scheme   nn_convfrq   = ', nn_convfrq
+         WRITE(numout,*) 'lim_hdf_init : Ice horizontal diffusion'
+         WRITE(numout,*) '~~~~~~~~~~~'
+         WRITE(numout,*) '   horizontal diffusivity calculation                          nn_ahi0      = ', nn_ahi0
+         WRITE(numout,*) '   horizontal diffusivity coeff. (orca2 grid)                  rn_ahi0_ref  = ', rn_ahi0_ref
+         WRITE(numout,*) '   convergence check frequency of the Crant-Nicholson scheme   nn_convfrq   = ', nn_convfrq
       ENDIF
+      !
+      !  Diffusion coefficients
+      SELECT CASE( nn_ahi0 )
+
+      CASE( -1 )
+         ahiu(:,:) = 0._wp
+         ahiv(:,:) = 0._wp
+
+         IF(lwp) WRITE(numout,*) ''
+         IF(lwp) WRITE(numout,*) '   No sea-ice diffusion applied'
+
+      CASE( 0 )
+         ahiu(:,:) = rn_ahi0_ref
+         ahiv(:,:) = rn_ahi0_ref
+
+         IF(lwp) WRITE(numout,*) ''
+         IF(lwp) WRITE(numout,*) '   laplacian operator: ahim constant = rn_ahi0_ref'
+
+      CASE( 1 ) 
+
+         zd_max = MAX( MAXVAL( e1t(:,:) ), MAXVAL( e2t(:,:) ) )
+         IF( lk_mpp )   CALL mpp_max( zd_max )          ! max over the global domain
+         
+         ahiu(:,:) = rn_ahi0_ref * zd_max * 1.e-05_wp   ! 1.e05 = 100km = max grid space at 60deg latitude in orca2
+                                                        !                    (60deg = min latitude for ice cover)  
+         ahiv(:,:) = rn_ahi0_ref * zd_max * 1.e-05_wp
+
+         IF(lwp) WRITE(numout,*) ''
+         IF(lwp) WRITE(numout,*) '   laplacian operator: ahim proportional to max of e1 e2 over the domain (', zd_max, ')'
+         IF(lwp) WRITE(numout,*) '   value for ahim = ', rn_ahi0_ref * zd_max * 1.e-05_wp 
+         
+      CASE( 2 ) 
+
+         zd_max = MAX( MAXVAL( e1t(:,:) ), MAXVAL( e2t(:,:) ) )
+         IF( lk_mpp )   CALL mpp_max( zd_max )   ! max over the global domain
+         
+         za00 = rn_ahi0_ref * 1.e-05_wp          ! 1.e05 = 100km = max grid space at 60deg latitude in orca2
+                                                 !                    (60deg = min latitude for ice cover)  
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               ahiu(ji,jj) = za00 * MAX( e1t(ji,jj), e2t(ji,jj) ) * umask(ji,jj,1)
+               ahiv(ji,jj) = za00 * MAX( e1f(ji,jj), e2f(ji,jj) ) * vmask(ji,jj,1)
+            END DO
+         END DO
+         !
+         IF(lwp) WRITE(numout,*) ''
+         IF(lwp) WRITE(numout,*) '   laplacian operator: ahim proportional to e1'
+         IF(lwp) WRITE(numout,*) '   maximum grid-spacing = ', zd_max, ' maximum value for ahim = ', za00*zd_max
+         
+      END SELECT
       !
    END SUBROUTINE lim_hdf_init
Index: branches/2015/nemo_v3_6_STABLE/NEMOGCM/NEMO/LIM_SRC_3/limthd.F90
===================================================================
--- branches/2015/nemo_v3_6_STABLE/NEMOGCM/NEMO/LIM_SRC_3/limthd.F90	(revision 7564)
+++ branches/2015/nemo_v3_6_STABLE/NEMOGCM/NEMO/LIM_SRC_3/limthd.F90	(revision 7607)
@@ -613,5 +613,4 @@
          CALL tab_1d_2d( nbpb, qns_ice(:,:,jl), npb, qns_ice_1d(1:nbpb) , jpi, jpj)
          CALL tab_1d_2d( nbpb, ftr_ice(:,:,jl), npb, ftr_ice_1d(1:nbpb) , jpi, jpj )
-         !         
       END SELECT
 
@@ -633,6 +632,6 @@
       INTEGER  ::   ios                 ! Local integer output status for namelist read
       NAMELIST/namicethd/ rn_hnewice, ln_frazil, rn_maxfrazb, rn_vfrazb, rn_Cfrazb,                       &
-         &                rn_himin, rn_betas, rn_kappa_i, nn_conv_dif, rn_terr_dif, nn_ice_thcon, &
-         &                nn_monocat, ln_it_qnsice
+         &                rn_himin, rn_betas, rn_kappa_i, nn_conv_dif, rn_terr_dif, nn_ice_thcon,         &
+         &                rn_cdsn, nn_monocat, ln_it_qnsice
       !!-------------------------------------------------------------------
       !
@@ -673,4 +672,5 @@
          WRITE(numout,*)'      maximal err. on T for heat diffusion computation        rn_terr_dif  = ', rn_terr_dif
          WRITE(numout,*)'      switch for comp. of thermal conductivity in the ice     nn_ice_thcon = ', nn_ice_thcon
+         WRITE(numout,*)'      thermal conductivity of the snow                        rn_cdsn      = ', rn_cdsn
          WRITE(numout,*)'      check heat conservation in the ice/snow                 con_i        = ', con_i
          WRITE(numout,*)'      virtual ITD mono-category parameterizations (1) or not  nn_monocat   = ', nn_monocat
Index: branches/2015/nemo_v3_6_STABLE/NEMOGCM/NEMO/LIM_SRC_3/limthd_dif.F90
===================================================================
--- branches/2015/nemo_v3_6_STABLE/NEMOGCM/NEMO/LIM_SRC_3/limthd_dif.F90	(revision 7564)
+++ branches/2015/nemo_v3_6_STABLE/NEMOGCM/NEMO/LIM_SRC_3/limthd_dif.F90	(revision 7607)
@@ -376,6 +376,6 @@
 
                ! Effective thickness he (zhe)
-               zfac     = 1._wp / ( rcdsn + zkimean )
-               zratio_s = rcdsn   * zfac
+               zfac     = 1._wp / ( rn_cdsn + zkimean )
+               zratio_s = rn_cdsn   * zfac
                zratio_i = zkimean * zfac
                zhe      = zratio_s * ht_i_1d(ji) + zratio_i * ht_s_1d(ji)
@@ -400,11 +400,11 @@
          DO ji = kideb, kiut
             zfac                  =  1. / MAX( epsi10 , zh_s(ji) )
-            zkappa_s(ji,0)        = zghe(ji) * rcdsn * zfac
-            zkappa_s(ji,nlay_s)   = zghe(ji) * rcdsn * zfac
+            zkappa_s(ji,0)        = zghe(ji) * rn_cdsn * zfac
+            zkappa_s(ji,nlay_s)   = zghe(ji) * rn_cdsn * zfac
          END DO
 
          DO jk = 1, nlay_s-1
             DO ji = kideb , kiut
-               zkappa_s(ji,jk)    = zghe(ji) * 2.0 * rcdsn / MAX( epsi10, 2.0 * zh_s(ji) )
+               zkappa_s(ji,jk)    = zghe(ji) * 2.0 * rn_cdsn / MAX( epsi10, 2.0 * zh_s(ji) )
             END DO
          END DO
@@ -422,6 +422,6 @@
             zkappa_i(ji,0)        = zghe(ji) * ztcond_i(ji,0) * zfac
             zkappa_i(ji,nlay_i)   = zghe(ji) * ztcond_i(ji,nlay_i) * zfac
-            zkappa_s(ji,nlay_s)   = zghe(ji) * zghe(ji) * 2.0 * rcdsn * ztcond_i(ji,0) / & 
-           &                        MAX( epsi10, ( zghe(ji) * ztcond_i(ji,0) * zh_s(ji) + zghe(ji) * rcdsn * zh_i(ji) ) )
+            zkappa_s(ji,nlay_s)   = zghe(ji) * zghe(ji) * 2.0 * rn_cdsn * ztcond_i(ji,0) / & 
+           &                        MAX( epsi10, ( zghe(ji) * ztcond_i(ji,0) * zh_s(ji) + zghe(ji) * rn_cdsn * zh_i(ji) ) )
             zkappa_i(ji,0)        = zkappa_s(ji,nlay_s) * isnow(ji) + zkappa_i(ji,0) * ( 1._wp - isnow(ji) )
          END DO
@@ -697,4 +697,5 @@
                &             ( isnow(ji) * t_s_1d(ji,1) + ( 1._wp - isnow(ji) ) * t_i_1d(ji,1) ) ) / zdiagbis(ji,numeqmin(ji))  
          END DO
+
          !
          !--------------------------------------------------------------------------
Index: branches/2015/nemo_v3_6_STABLE/NEMOGCM/NEMO/OPA_SRC/DOM/phycst.F90
===================================================================
--- branches/2015/nemo_v3_6_STABLE/NEMOGCM/NEMO/OPA_SRC/DOM/phycst.F90	(revision 7564)
+++ branches/2015/nemo_v3_6_STABLE/NEMOGCM/NEMO/OPA_SRC/DOM/phycst.F90	(revision 7607)
@@ -65,7 +65,6 @@
 #if defined key_lim3 || defined key_cice
    REAL(wp), PUBLIC ::   rhoic    =  917._wp         !: volumic mass of sea ice                               [kg/m3]
-   REAL(wp), PUBLIC ::   rcdic    =    2.034396_wp   !: thermal conductivity of fresh ice
-   REAL(wp), PUBLIC ::   rcdsn    =    0.31_wp       !: thermal conductivity of snow
-   REAL(wp), PUBLIC ::   cpic     = 2067.0_wp        !: specific heat for ice 
+   REAL(wp), PUBLIC ::   rcdic    =    2.034396_wp   !: thermal conductivity of fresh ice                     [W/m/K]
+   REAL(wp), PUBLIC ::   cpic     = 2067.0_wp        !: specific heat of fresh ice                            [J/kg/K]
    REAL(wp), PUBLIC ::   lsub     =    2.834e+6_wp   !: pure ice latent heat of sublimation                   [J/kg]
    REAL(wp), PUBLIC ::   lfus     =    0.334e+6_wp   !: latent heat of fusion of fresh ice                    [J/kg]
@@ -83,4 +82,7 @@
    REAL(wp), PUBLIC ::   xlic     =  300.33e+6_wp    !: volumetric latent heat fusion of ice                  [J/m3]
    REAL(wp), PUBLIC ::   xsn      =    2.8e+6_wp     !: volumetric latent heat of sublimation of snow         [J/m3]
+#endif
+#if defined key_cice
+   REAL(wp), PUBLIC ::   rcdsn    =    0.31_wp       !: thermal conductivity of snow                          [W/m/K], now namelist parameter for LIM3
 #endif
 #if defined key_lim3
@@ -177,6 +179,8 @@
       IF(lwp) THEN
          WRITE(numout,*)
+#if defined key_cice
          WRITE(numout,*) '          thermal conductivity of the snow          = ', rcdsn   , ' J/s/m/K'
-         WRITE(numout,*) '          thermal conductivity of the ice           = ', rcdic   , ' J/s/m/K'
+#endif
+         WRITE(numout,*) '          thermal conductivity of pure ice          = ', rcdic   , ' J/s/m/K'
          WRITE(numout,*) '          fresh ice specific heat                   = ', cpic    , ' J/kg/K'
          WRITE(numout,*) '          latent heat of fusion of fresh ice / snow = ', lfus    , ' J/kg'
Index: branches/2015/nemo_v3_6_STABLE/NEMOGCM/NEMO/OPA_SRC/SBC/albedo.F90
===================================================================
--- branches/2015/nemo_v3_6_STABLE/NEMOGCM/NEMO/OPA_SRC/SBC/albedo.F90	(revision 7564)
+++ branches/2015/nemo_v3_6_STABLE/NEMOGCM/NEMO/OPA_SRC/SBC/albedo.F90	(revision 7607)
@@ -39,5 +39,5 @@
    !                             !!* namelist namsbc_alb
    INTEGER  ::   nn_ice_alb
-   REAL(wp) ::   rn_albice
+   REAL(wp) ::   rn_alb_sdry, rn_alb_smlt, rn_alb_idry, rn_alb_imlt
 
    !!----------------------------------------------------------------------
@@ -101,4 +101,8 @@
       IF( albd_init == 0 )   CALL albedo_init      ! initialization 
 
+      ralb_sf = rn_alb_sdry ! dry snow
+      ralb_sm = rn_alb_smlt ! melting snow
+      ralb_if = rn_alb_idry ! bare frozen ice
+      ralb_im = rn_alb_imlt ! bare puddled ice 
       
       SELECT CASE ( nn_ice_alb )
@@ -109,9 +113,9 @@
       CASE( 0 )
        
-         ralb_sf = 0.80       ! dry snow
-         ralb_sm = 0.65       ! melting snow
-         ralb_if = 0.72       ! bare frozen ice
-         ralb_im = rn_albice  ! bare puddled ice 
-         
+         !ralb_sf = 0.80       ! dry snow
+         !ralb_sm = 0.65       ! melting snow
+         !ralb_if = 0.72       ! bare frozen ice
+         !ralb_im = ...        ! bare puddled ice 
+
          !  Computation of ice albedo (free of snow)
          WHERE     ( ph_snw == 0._wp .AND. pt_ice >= rt0_ice )   ;   zalb(:,:,:) = ralb_im
@@ -163,9 +167,9 @@
       CASE( 1 ) 
 
-         ralb_im = rn_albice  ! bare puddled ice
+!        ralb_im = ...        ! bare puddled ice
 ! compilation of values from literature
-         ralb_sf = 0.85      ! dry snow
-         ralb_sm = 0.75      ! melting snow
-         ralb_if = 0.60      ! bare frozen ice
+!        ralb_sf = 0.85      ! dry snow
+!        ralb_sm = 0.75      ! melting snow
+!        ralb_if = 0.60      ! bare frozen ice
 ! Perovich et al 2002 (Sheba) => the only dataset for which all types of ice/snow were retrieved
 !         ralb_sf = 0.85       ! dry snow
@@ -248,5 +252,5 @@
       !!----------------------------------------------------------------------
       INTEGER  ::   ios                 ! Local integer output status for namelist read
-      NAMELIST/namsbc_alb/ nn_ice_alb, rn_albice 
+      NAMELIST/namsbc_alb/ nn_ice_alb, rn_alb_sdry, rn_alb_smlt, rn_alb_idry , rn_alb_imlt
       !!----------------------------------------------------------------------
       !
@@ -268,5 +272,8 @@
          WRITE(numout,*) '   Namelist namsbc_alb : albedo '
          WRITE(numout,*) '      choose the albedo parameterization                  nn_ice_alb = ', nn_ice_alb
-         WRITE(numout,*) '      albedo of bare puddled ice                          rn_albice  = ', rn_albice
+         WRITE(numout,*) '      albedo of dry snow                                  rn_alb_sdry = ', rn_alb_sdry
+         WRITE(numout,*) '      albedo of melting snow                              rn_alb_smlt = ', rn_alb_smlt
+         WRITE(numout,*) '      albedo of dry ice                                   rn_alb_idry = ', rn_alb_idry
+         WRITE(numout,*) '      albedo of bare puddled ice                          rn_alb_imlt = ', rn_alb_imlt
       ENDIF
       !
Index: branches/2015/nemo_v3_6_STABLE/NEMOGCM/NEMO/OPA_SRC/SBC/sbc_oce.F90
===================================================================
--- branches/2015/nemo_v3_6_STABLE/NEMOGCM/NEMO/OPA_SRC/SBC/sbc_oce.F90	(revision 7564)
+++ branches/2015/nemo_v3_6_STABLE/NEMOGCM/NEMO/OPA_SRC/SBC/sbc_oce.F90	(revision 7607)
@@ -113,4 +113,5 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   rnf    , rnf_b    !: river runoff        [Kg/m2/s]  
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   fwfisf , fwfisf_b !: ice shelf melting   [Kg/m2/s]  
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   fwficb , fwficb_b !: iceberg melting [Kg/m2/s]  
    !!
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::  sbc_tsc, sbc_tsc_b  !: sbc content trend                      [K.m/s] jpi,jpj,jpts
@@ -164,5 +165,6 @@
          !
       ALLOCATE( fwfisf  (jpi,jpj), rnf  (jpi,jpj) , sbc_tsc  (jpi,jpj,jpts) , qsr_hc  (jpi,jpj,jpk) ,     &
-         &      fwfisf_b(jpi,jpj), rnf_b(jpi,jpj) , sbc_tsc_b(jpi,jpj,jpts) , qsr_hc_b(jpi,jpj,jpk) , STAT=ierr(3) )
+         &      fwfisf_b(jpi,jpj), rnf_b(jpi,jpj) , sbc_tsc_b(jpi,jpj,jpts) , qsr_hc_b(jpi,jpj,jpk) ,     &
+         &      fwficb  (jpi,jpj), fwficb_b(jpi,jpj), STAT=ierr(3) )
          !
       ALLOCATE( tprecip(jpi,jpj) , sprecip(jpi,jpj) , fr_i(jpi,jpj) ,     &
Index: branches/2015/nemo_v3_6_STABLE/NEMOGCM/NEMO/OPA_SRC/SBC/sbccpl.F90
===================================================================
--- branches/2015/nemo_v3_6_STABLE/NEMOGCM/NEMO/OPA_SRC/SBC/sbccpl.F90	(revision 7564)
+++ branches/2015/nemo_v3_6_STABLE/NEMOGCM/NEMO/OPA_SRC/SBC/sbccpl.F90	(revision 7607)
@@ -43,4 +43,5 @@
    USE eosbn2
    USE sbcrnf   , ONLY : l_rnfcpl
+   USE sbcisf   , ONLY : l_isfcpl
 #if defined key_cpl_carbon_cycle
    USE p4zflx, ONLY : oce_co2
@@ -105,5 +106,7 @@
    INTEGER, PARAMETER ::   jpr_e3t1st = 41            ! first T level thickness 
    INTEGER, PARAMETER ::   jpr_fraqsr = 42            ! fraction of solar net radiation absorbed in the first ocean level
-   INTEGER, PARAMETER ::   jprcv      = 42            ! total number of fields received
+   INTEGER, PARAMETER ::   jpr_isf    = 43
+   INTEGER, PARAMETER ::   jpr_icb    = 44
+   INTEGER, PARAMETER ::   jprcv      = 44            ! total number of fields received
 
    INTEGER, PARAMETER ::   jps_fice   =  1            ! ice fraction sent to the atmosphere
@@ -149,5 +152,5 @@
    ! Received from the atmosphere                     !
    TYPE(FLD_C) ::   sn_rcv_w10m, sn_rcv_taumod, sn_rcv_tau, sn_rcv_dqnsdt, sn_rcv_qsr, sn_rcv_qns, sn_rcv_emp, sn_rcv_rnf
-   TYPE(FLD_C) ::   sn_rcv_cal, sn_rcv_iceflx, sn_rcv_co2                        
+   TYPE(FLD_C) ::   sn_rcv_cal, sn_rcv_iceflx, sn_rcv_co2, sn_rcv_icb, sn_rcv_isf                               
    ! Other namelist parameters                        !
    INTEGER     ::   nn_cplmodel            ! Maximum number of models to/from which NEMO is potentialy sending/receiving data
@@ -219,5 +222,5 @@
          &                  sn_rcv_w10m, sn_rcv_taumod, sn_rcv_tau  , sn_rcv_dqnsdt, sn_rcv_qsr,      &
          &                  sn_rcv_qns , sn_rcv_emp   , sn_rcv_rnf  , sn_rcv_cal   , sn_rcv_iceflx,   &
-         &                  sn_rcv_co2 , nn_cplmodel  , ln_usecplmask
+         &                  sn_rcv_co2 , sn_rcv_icb , sn_rcv_isf, nn_cplmodel  , ln_usecplmask
       !!---------------------------------------------------------------------
       !
@@ -258,4 +261,6 @@
          WRITE(numout,*)'      runoffs                         = ', TRIM(sn_rcv_rnf%cldes   ), ' (', TRIM(sn_rcv_rnf%clcat   ), ')'
          WRITE(numout,*)'      calving                         = ', TRIM(sn_rcv_cal%cldes   ), ' (', TRIM(sn_rcv_cal%clcat   ), ')'
+         WRITE(numout,*)'      iceberg                         = ', TRIM(sn_rcv_icb%cldes   ), ' (', TRIM(sn_rcv_icb%clcat   ), ')'
+         WRITE(numout,*)'      ice shelf                       = ', TRIM(sn_rcv_isf%cldes   ), ' (', TRIM(sn_rcv_isf%clcat   ), ')'
          WRITE(numout,*)'      sea ice heat fluxes             = ', TRIM(sn_rcv_iceflx%cldes), ' (', TRIM(sn_rcv_iceflx%clcat), ')'
          WRITE(numout,*)'      atm co2                         = ', TRIM(sn_rcv_co2%cldes   ), ' (', TRIM(sn_rcv_co2%clcat   ), ')'
@@ -397,7 +402,8 @@
       END SELECT
 
-      !                                                      ! ------------------------- !
-      !                                                      !     Runoffs & Calving     !   
-      !                                                      ! ------------------------- !
+
+      !                                                      ! ---------------------------------------------------- !
+      !                                                      !     Runoffs, Calving, Iceberg, Iceshelf cavities     !   
+      !                                                      ! ---------------------------------------------------- !
       srcv(jpr_rnf   )%clname = 'O_Runoff'
       IF( TRIM( sn_rcv_rnf%cldes ) == 'coupled' ) THEN
@@ -409,5 +415,13 @@
       ENDIF
       !
-      srcv(jpr_cal   )%clname = 'OCalving'   ;   IF( TRIM( sn_rcv_cal%cldes ) == 'coupled' )   srcv(jpr_cal)%laction = .TRUE.
+      srcv(jpr_cal)%clname = 'OCalving'   ;  IF( TRIM( sn_rcv_cal%cldes) == 'coupled' )   srcv(jpr_cal)%laction = .TRUE.
+      srcv(jpr_isf)%clname = 'OIcshelf'   ;  IF( TRIM( sn_rcv_isf%cldes) == 'coupled' )   srcv(jpr_isf)%laction = .TRUE.
+      srcv(jpr_icb)%clname = 'OIceberg'   ;  IF( TRIM( sn_rcv_icb%cldes) == 'coupled' )   srcv(jpr_icb)%laction = .TRUE.
+
+      IF( srcv(jpr_isf)%laction .AND. nn_isf > 0 ) THEN
+         l_isfcpl             = .TRUE.                      ! -> no need to read isf in sbcisf
+         IF(lwp) WRITE(numout,*)
+         IF(lwp) WRITE(numout,*) '   iceshelf received from oasis '
+      ENDIF
 
       !                                                      ! ------------------------- !
@@ -1071,7 +1085,14 @@
          ENDIF
          !
+         !   
          !                                                        ! runoffs and calving (added in emp)
-         IF( srcv(jpr_rnf)%laction )     rnf(:,:) = frcv(jpr_rnf)%z3(:,:,1)
+         IF( srcv(jpr_rnf)%laction )     rnf(:,:)  = frcv(jpr_rnf)%z3(:,:,1)
          IF( srcv(jpr_cal)%laction )     zemp(:,:) = zemp(:,:) - frcv(jpr_cal)%z3(:,:,1)
+
+         IF( srcv(jpr_icb)%laction )  THEN 
+             fwficb(:,:) = frcv(jpr_icb)%z3(:,:,1)
+             rnf(:,:)    = rnf(:,:) + fwficb(:,:)   ! iceberg added to runfofs
+         ENDIF
+         IF( srcv(jpr_isf)%laction )  fwfisf(:,:) = - frcv(jpr_isf)%z3(:,:,1)  ! fresh water flux from the isf (fwfisf <0 mean melting)  
          
          IF( ln_mixcpl ) THEN   ;   emp(:,:) = emp(:,:) * xcplmask(:,:,0) + zemp(:,:) * zmsk(:,:)
@@ -1091,4 +1112,7 @@
             ENDIF
          ENDIF
+         !
+         IF( srcv(jpr_icb)%laction )  zqns(:,:) = zqns(:,:) - frcv(jpr_icb)%z3(:,:,1) * lfus ! remove heat content associated to iceberg melting
+         !
          IF( ln_mixcpl ) THEN   ;   qns(:,:) = qns(:,:) * xcplmask(:,:,0) + zqns(:,:) * zmsk(:,:)
          ELSE                   ;   qns(:,:) =                              zqns(:,:)
@@ -1468,4 +1492,15 @@
       ENDIF
 
+      IF( srcv(jpr_icb)%laction )  THEN 
+         fwficb(:,:) = frcv(jpr_icb)%z3(:,:,1)
+         rnf(:,:)    = rnf(:,:) + fwficb(:,:)   ! iceberg added to runoffs
+         CALL iom_put( 'iceberg_cea', frcv(jpr_icb)%z3(:,:,1) )
+      ENDIF
+      IF( srcv(jpr_isf)%laction )  THEN
+        fwfisf(:,:) = - frcv(jpr_isf)%z3(:,:,1)  ! fresh water flux from the isf (fwfisf <0 mean melting)  
+        CALL iom_put( 'iceshelf_cea', frcv(jpr_isf)%z3(:,:,1) )
+      ENDIF
+
+
       IF( ln_mixcpl ) THEN
          emp_tot(:,:) = emp_tot(:,:) * xcplmask(:,:,0) + zemp_tot(:,:) * zmsk(:,:)
@@ -1504,4 +1539,16 @@
       ENDIF
 
+
+      IF( srcv(jpr_icb)%laction )  THEN 
+         fwficb(:,:) = frcv(jpr_icb)%z3(:,:,1)
+         rnf(:,:)    = rnf(:,:) + fwficb(:,:)   ! iceberg added to runoffs
+         CALL iom_put( 'iceberg_cea', frcv(jpr_icb)%z3(:,:,1) )
+      ENDIF
+      IF( srcv(jpr_isf)%laction )  THEN
+        fwfisf(:,:) = - frcv(jpr_isf)%z3(:,:,1)  ! fresh water flux from the isf (fwfisf <0 mean melting)  
+        CALL iom_put( 'iceshelf_cea', frcv(jpr_isf)%z3(:,:,1) )
+      ENDIF
+
+
       IF( ln_mixcpl ) THEN
          emp_tot(:,:) = emp_tot(:,:) * xcplmask(:,:,0) + zemp_tot(:,:) * zmsk(:,:)
@@ -1570,4 +1617,11 @@
          IF( iom_use('hflx_cal_cea') )   CALL iom_put( 'hflx_cal_cea', - frcv(jpr_cal)%z3(:,:,1) * lfus )   ! heat flux from calving
       ENDIF
+
+!!chris     
+!!    The heat content associated to the ice shelf in removed in the routine sbcisf.F90
+      !
+      IF( srcv(jpr_icb)%laction )  zqns_tot(:,:) = zqns_tot(:,:) - frcv(jpr_icb)%z3(:,:,1) * lfus ! remove heat content associated to iceberg melting
+      !
+!!      !
 
 #if defined key_lim3      
Index: branches/2015/nemo_v3_6_STABLE/NEMOGCM/NEMO/OPA_SRC/SBC/sbcisf.F90
===================================================================
--- branches/2015/nemo_v3_6_STABLE/NEMOGCM/NEMO/OPA_SRC/SBC/sbcisf.F90	(revision 7564)
+++ branches/2015/nemo_v3_6_STABLE/NEMOGCM/NEMO/OPA_SRC/SBC/sbcisf.F90	(revision 7607)
@@ -55,4 +55,6 @@
    INTEGER,    PUBLIC, ALLOCATABLE, SAVE, DIMENSION (:,:)     ::  misfkt, misfkb         !:Level of ice shelf base
 
+   LOGICAL, PUBLIC ::   l_isfcpl = .false.       ! isf recieved from oasis
+
 
    REAL(wp), PUBLIC, SAVE ::   rcpi   = 2000.0_wp     ! phycst ?
@@ -94,4 +96,6 @@
     REAL(wp), DIMENSION(:,:,:), POINTER :: zfwfisf3d, zqhcisf3d, zqlatisf3d
     REAL(wp), DIMENSION(:,:  ), POINTER :: zqhcisf2d
+    REAL(wp)                            :: zhisf
+
       !
       !!---------------------------------------------------------------------
@@ -142,7 +146,9 @@
             misfkt(:,:)    = mikt(:,:)         ! same indice for bg03 et cav => used in isfdiv
          ELSE IF ((nn_isf == 3) .OR. (nn_isf == 2)) THEN
-            ALLOCATE( sf_rnfisf(1), STAT=ierror )
-            ALLOCATE( sf_rnfisf(1)%fnow(jpi,jpj,1), sf_rnfisf(1)%fdta(jpi,jpj,1,2) )
-            CALL fld_fill( sf_rnfisf, (/ sn_rnfisf /), cn_dirisf, 'sbc_isf_init', 'read fresh water flux isf data', 'namsbc_isf' )
+            IF( .NOT.l_isfcpl ) THEN
+               ALLOCATE( sf_rnfisf(1), STAT=ierror )
+               ALLOCATE( sf_rnfisf(1)%fnow(jpi,jpj,1), sf_rnfisf(1)%fdta(jpi,jpj,1,2) )
+               CALL fld_fill( sf_rnfisf, (/ sn_rnfisf /), cn_dirisf, 'sbc_isf_init', 'read fresh water flux isf data', 'namsbc_isf' )
+             ENDIF
 
             !: read effective lenght (BG03)
@@ -185,9 +191,11 @@
             
             ! load variable used in fldread (use for temporal interpolation of isf fwf forcing)
-            ALLOCATE( sf_fwfisf(1), sf_qisf(1), STAT=ierror )
-            ALLOCATE( sf_fwfisf(1)%fnow(jpi,jpj,1), sf_fwfisf(1)%fdta(jpi,jpj,1,2) )
-            ALLOCATE( sf_qisf(1)%fnow(jpi,jpj,1), sf_qisf(1)%fdta(jpi,jpj,1,2) )
-            CALL fld_fill( sf_fwfisf, (/ sn_fwfisf /), cn_dirisf, 'sbc_isf_init', 'read fresh water flux isf data', 'namsbc_isf' )
-            !CALL fld_fill( sf_qisf  , (/ sn_qisf   /), cn_dirisf, 'sbc_isf_init', 'read heat flux isf data'       , 'namsbc_isf' )
+            IF( .NOT.l_isfcpl ) THEN
+               ALLOCATE( sf_fwfisf(1), sf_qisf(1), STAT=ierror )
+               ALLOCATE( sf_fwfisf(1)%fnow(jpi,jpj,1), sf_fwfisf(1)%fdta(jpi,jpj,1,2) )
+               ALLOCATE( sf_qisf(1)%fnow(jpi,jpj,1), sf_qisf(1)%fdta(jpi,jpj,1,2) )
+               CALL fld_fill( sf_fwfisf, (/ sn_fwfisf /), cn_dirisf, 'sbc_isf_init', 'read fresh water flux isf data', 'namsbc_isf' )
+               !CALL fld_fill( sf_qisf  , (/ sn_qisf   /), cn_dirisf, 'sbc_isf_init', 'read heat flux isf data'       , 'namsbc_isf' )
+            ENDIF
          END IF
          
@@ -242,5 +250,7 @@
             CALL iom_put('vtbl',vtbl(:,:))
             ! compute fwf and heat flux
-            CALL sbc_isf_cav (kt)
+            IF( .NOT.l_isfcpl ) THEN    ;   CALL sbc_isf_cav (kt)
+            ELSE                        ;   qisf(:,:)  = fwfisf(:,:) * lfusisf              ! heat        flux
+            ENDIF
 
          ELSE IF (nn_isf == 2) THEN
@@ -251,6 +261,8 @@
          ELSE IF (nn_isf == 3) THEN
             ! specified runoff in depth (Mathiot et al., XXXX in preparation)
-            CALL fld_read ( kt, nn_fsbc, sf_rnfisf   )
-            fwfisf(:,:) = - sf_rnfisf(1)%fnow(:,:,1)         ! fresh water flux from the isf (fwfisf <0 mean melting) 
+            IF( .NOT.l_isfcpl ) THEN
+               CALL fld_read ( kt, nn_fsbc, sf_rnfisf   )
+               fwfisf(:,:) = - sf_rnfisf(1)%fnow(:,:,1)         ! fresh water flux from the isf (fwfisf <0 mean melting) 
+            ENDIF
             qisf(:,:)   = fwfisf(:,:) * lfusisf              ! heat        flux
             stbl(:,:)   = soce
@@ -258,7 +270,9 @@
          ELSE IF (nn_isf == 4) THEN
             ! specified fwf and heat flux forcing beneath the ice shelf
-            CALL fld_read ( kt, nn_fsbc, sf_fwfisf   )
-            !CALL fld_read ( kt, nn_fsbc, sf_qisf   )
-            fwfisf(:,:) = sf_fwfisf(1)%fnow(:,:,1)            ! fwf
+            IF( .NOT.l_isfcpl ) THEN
+               CALL fld_read ( kt, nn_fsbc, sf_fwfisf   )
+               !CALL fld_read ( kt, nn_fsbc, sf_qisf   )
+               fwfisf(:,:) = sf_fwfisf(1)%fnow(:,:,1)            ! fwf
+            ENDIF
             qisf(:,:)   = fwfisf(:,:) * lfusisf              ! heat        flux
             !qisf(:,:)   = sf_qisf(1)%fnow(:,:,1)              ! heat flux
@@ -288,14 +302,15 @@
          CALL lbc_lnk(qisf(:,:)     ,'T',1.)
 
-!=============================================================================================================================================
-         IF ( iom_use('fwfisf3d') .OR. iom_use('qlatisf3d') .OR. iom_use('qhcisf3d') .OR. iom_use('qhcisf')) THEN
+         ! Diagnostics
+         IF( iom_use('fwfisf3d') .OR. iom_use('qlatisf3d') .OR. iom_use('qhcisf3d') .OR. iom_use('qhcisf')) THEN
+            !
             CALL wrk_alloc( jpi,jpj,jpk, zfwfisf3d, zqhcisf3d, zqlatisf3d )
             CALL wrk_alloc( jpi,jpj,     zqhcisf2d                        )
-
+            !
             zfwfisf3d(:,:,:) = 0.0_wp                         ! 3d ice shelf melting (kg/m2/s)
             zqhcisf3d(:,:,:) = 0.0_wp                         ! 3d heat content flux (W/m2)
             zqlatisf3d(:,:,:)= 0.0_wp                         ! 3d ice shelf melting latent heat flux (W/m2)
             zqhcisf2d(:,:)   = fwfisf(:,:) * zt_frz * rcp     ! 2d heat content flux (W/m2)
-
+            !
             DO jj = 1,jpj
                DO ji = 1,jpi
@@ -303,24 +318,27 @@
                   ikb = misfkb(ji,jj)
                   DO jk = ikt, ikb - 1
-                     zfwfisf3d (ji,jj,jk) = zfwfisf3d (ji,jj,jk) + fwfisf   (ji,jj) * r1_hisf_tbl(ji,jj) * fse3t(ji,jj,jk)
-                     zqhcisf3d (ji,jj,jk) = zqhcisf3d (ji,jj,jk) + zqhcisf2d(ji,jj) * r1_hisf_tbl(ji,jj) * fse3t(ji,jj,jk)
-                     zqlatisf3d(ji,jj,jk) = zqlatisf3d(ji,jj,jk) + qisf     (ji,jj) * r1_hisf_tbl(ji,jj) * fse3t(ji,jj,jk)
+                     zhisf = r1_hisf_tbl(ji,jj) * fse3t(ji,jj,jk)
+                     zfwfisf3d (ji,jj,jk) = zfwfisf3d (ji,jj,jk) + fwfisf(ji,jj)    * zhisf
+                     zqhcisf3d (ji,jj,jk) = zqhcisf3d (ji,jj,jk) + zqhcisf2d(ji,jj) * zhisf
+                     zqlatisf3d(ji,jj,jk) = zqlatisf3d(ji,jj,jk) + qisf(ji,jj)      * zhisf
                   END DO
-                  zfwfisf3d (ji,jj,jk) = zfwfisf3d (ji,jj,jk) + fwfisf   (ji,jj) * r1_hisf_tbl(ji,jj) * ralpha(ji,jj) * fse3t(ji,jj,jk)
-                  zqhcisf3d (ji,jj,jk) = zqhcisf3d (ji,jj,jk) + zqhcisf2d(ji,jj) * r1_hisf_tbl(ji,jj) * ralpha(ji,jj) * fse3t(ji,jj,jk)
-                  zqlatisf3d(ji,jj,jk) = zqlatisf3d(ji,jj,jk) + qisf     (ji,jj) * r1_hisf_tbl(ji,jj) * ralpha(ji,jj) * fse3t(ji,jj,jk)
+                  jk = ikb
+                  zhisf = r1_hisf_tbl(ji,jj) * fse3t(ji,jj,jk)
+                  zfwfisf3d (ji,jj,jk) = zfwfisf3d (ji,jj,jk) + fwfisf   (ji,jj) * zhisf * ralpha(ji,jj) 
+                  zqhcisf3d (ji,jj,jk) = zqhcisf3d (ji,jj,jk) + zqhcisf2d(ji,jj) * zhisf * ralpha(ji,jj)
+                  zqlatisf3d(ji,jj,jk) = zqlatisf3d(ji,jj,jk) + qisf     (ji,jj) * zhisf * ralpha(ji,jj)
                END DO
             END DO
-
-            CALL iom_put('fwfisf3d' , zfwfisf3d (:,:,:))
-            CALL iom_put('qlatisf3d', zqlatisf3d(:,:,:))
-            CALL iom_put('qhcisf3d' , zqhcisf3d (:,:,:))
-            CALL iom_put('qhcisf'   , zqhcisf2d (:,:  ))
-
+            !
+            CALL iom_put( 'fwfisf3d' , zfwfisf3d (:,:,:) )
+            CALL iom_put( 'qlatisf3d', zqlatisf3d(:,:,:) )
+            CALL iom_put( 'qhcisf3d' , zqhcisf3d (:,:,:) )
+            CALL iom_put( 'qhcisf'   , zqhcisf2d (:,:  ) )
+            !
             CALL wrk_dealloc( jpi,jpj,jpk, zfwfisf3d, zqhcisf3d, zqlatisf3d )
             CALL wrk_dealloc( jpi,jpj,     zqhcisf2d                        )
+            !
          END IF
-!=============================================================================================================================================
-
+         !
          IF( kt == nit000 ) THEN                          !   set the forcing field at nit000 - 1    !
             IF( ln_rstart .AND.    &                     ! Restart: read in restart file
Index: branches/2015/nemo_v3_6_STABLE/NEMOGCM/NEMO/OPA_SRC/stpctl.F90
===================================================================
--- branches/2015/nemo_v3_6_STABLE/NEMOGCM/NEMO/OPA_SRC/stpctl.F90	(revision 7564)
+++ branches/2015/nemo_v3_6_STABLE/NEMOGCM/NEMO/OPA_SRC/stpctl.F90	(revision 7607)
@@ -180,6 +180,6 @@
       ENDIF
 
-9200  FORMAT('it:', i8, ' iter:', i4, ' r: ',e16.10, ' b: ',e16.10 )
-9300  FORMAT(' it :', i8, ' ssh2: ', e16.10, ' Umax: ',e16.10,' Smin: ',e16.10)
+9200  FORMAT('it:', i8, ' iter:', i4, ' r: ',d23.16, ' b: ',d23.16 )
+9300  FORMAT(' it :', i8, ' ssh2: ', d23.16, ' Umax: ',d23.16,' Smin: ',d23.16)
       !
    END SUBROUTINE stp_ctl
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 7564)
+++ 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 7564)
+++ 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 7564)
+++ 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 7564)
+++ 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 7564)
+++ 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 7564)
+++ 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
