Index: branches/2016/dev_CNRS_2016/NEMOGCM/CONFIG/SHARED/namelist_pisces_ref
===================================================================
--- branches/2016/dev_CNRS_2016/NEMOGCM/CONFIG/SHARED/namelist_pisces_ref	(revision 7280)
+++ branches/2016/dev_CNRS_2016/NEMOGCM/CONFIG/SHARED/namelist_pisces_ref	(revision 7282)
@@ -66,4 +66,5 @@
    qdfelim    =  7.E-6    ! Optimal quota of diatoms
    caco3r     =  0.3      ! mean rain ratio
+   oxymin    =  1.E-6     ! Half-saturation constant for anoxia
 /
 !'''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
@@ -162,5 +163,4 @@
    xsiremlab =  0.03      ! fast remineralization rate of Si
    xsilab    =  0.5       ! Fraction of labile biogenic silica
-   oxymin    =  1.E-6     ! Half-saturation constant for anoxia
 /
 !'''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
Index: branches/2016/dev_CNRS_2016/NEMOGCM/NEMO/NST_SRC/agrif_lim2_interp.F90
===================================================================
--- branches/2016/dev_CNRS_2016/NEMOGCM/NEMO/NST_SRC/agrif_lim2_interp.F90	(revision 7280)
+++ branches/2016/dev_CNRS_2016/NEMOGCM/NEMO/NST_SRC/agrif_lim2_interp.F90	(revision 7282)
@@ -392,5 +392,4 @@
       INTEGER :: ji,jj,jn
       REAL(wp) :: zalpha
-      REAL(wp), DIMENSION(jpi,jpj,7) :: tabice_agr 
       !!-----------------------------------------------------------------------      
       !
@@ -529,4 +528,10 @@
             END DO
          END DO
+      ELSE
+         DO jj=MAX(j1,2),j2
+            DO ji=MAX(i1,2),i2
+               uice_agr(ji,jj) = tabres(ji,jj)
+            END DO
+         END DO
       ENDIF
 #else
@@ -541,4 +546,10 @@
             END DO
          END DO
+      ELSE
+         DO jj= j1, j2
+            DO ji= i1, i2
+               uice_agr(ji,jj) = tabres(ji,jj)
+            END DO
+         END DO
       ENDIF
 #endif
@@ -566,4 +577,10 @@
                   tabres(ji,jj) = e1f(ji-1,jj-1) * v_ice(ji,jj)
                ENDIF
+            END DO
+         END DO
+      ELSE
+         DO jj=MAX(j1,2),j2
+            DO ji=MAX(i1,2),i2
+               vice_agr(ji,jj) = tabres(ji,jj)
             END DO
          END DO
@@ -580,4 +597,10 @@
 	         END DO
 	      END DO
+      ELSE
+	      DO jj= j1 ,j2
+	         DO ji = i1, i2
+	            vice_agr(ji,jj) = tabres(ji,jj)
+	         END DO
+	      END DO
       ENDIF
 #endif
@@ -585,5 +608,5 @@
 
 
-   SUBROUTINE interp_adv_ice( tabres, i1, i2, j1, j2, before )
+   SUBROUTINE interp_adv_ice( tabres, i1, i2, j1, j2, k1, k2, before )
       !!-----------------------------------------------------------------------
       !!                    *** ROUTINE interp_adv_ice ***                           
@@ -593,6 +616,6 @@
       !!              put -9999 where no ice for correct extrapolation             
       !!-----------------------------------------------------------------------
-      INTEGER, INTENT(in) :: i1, i2, j1, j2
-      REAL(wp), DIMENSION(i1:i2,j1:j2,7), INTENT(inout) :: tabres
+      INTEGER, INTENT(in) :: i1, i2, j1, j2, k1, k2
+      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2), INTENT(inout) :: tabres
       LOGICAL, INTENT(in) :: before
       !!
@@ -601,19 +624,27 @@
       !
       IF( before ) THEN
-	      DO jj=j1,j2
-	         DO ji=i1,i2
-	            IF( tms(ji,jj) == 0. ) THEN
-	               tabres(ji,jj,:) = -9999. 
-	            ELSE
-	               tabres(ji,jj, 1) = frld  (ji,jj)
-	               tabres(ji,jj, 2) = hicif (ji,jj)
-	               tabres(ji,jj, 3) = hsnif (ji,jj)
-	               tabres(ji,jj, 4) = tbif  (ji,jj,1)
-	               tabres(ji,jj, 5) = tbif  (ji,jj,2)
-	               tabres(ji,jj, 6) = tbif  (ji,jj,3)
-	               tabres(ji,jj, 7) = qstoif(ji,jj)
-	            ENDIF
-	         END DO
-	      END DO
+         DO jj=j1,j2
+	    DO ji=i1,i2
+	       IF( tms(ji,jj) == 0. ) THEN
+	          tabres(ji,jj,:) = -9999 
+	       ELSE
+	          tabres(ji,jj, 1) = frld  (ji,jj)
+	          tabres(ji,jj, 2) = hicif (ji,jj)
+	          tabres(ji,jj, 3) = hsnif (ji,jj)
+	          tabres(ji,jj, 4) = tbif  (ji,jj,1)
+	          tabres(ji,jj, 5) = tbif  (ji,jj,2)
+	          tabres(ji,jj, 6) = tbif  (ji,jj,3)
+	          tabres(ji,jj, 7) = qstoif(ji,jj)
+	       ENDIF
+	    END DO
+         END DO
+      ELSE
+	 DO jj=j1,j2
+	    DO ji=i1,i2
+               DO jk=k1, k2
+	          tabice_agr(ji,jj,jk) = tabres(ji,jj,jk)
+               END DO
+	    END DO
+	 END DO
       ENDIF
       !
Index: branches/2016/dev_CNRS_2016/NEMOGCM/NEMO/OPA_SRC/DOM/domain.F90
===================================================================
--- branches/2016/dev_CNRS_2016/NEMOGCM/NEMO/OPA_SRC/DOM/domain.F90	(revision 7280)
+++ branches/2016/dev_CNRS_2016/NEMOGCM/NEMO/OPA_SRC/DOM/domain.F90	(revision 7282)
@@ -153,15 +153,15 @@
       !
          !       before        !          now          !       after         !
-         ;  gdept_b = gdept_0  ;   gdept_n = gdept_0   !        ---          ! depth of grid-points
-         ;  gdepw_b = gdepw_0  ;   gdepw_n = gdepw_0   !        ---          !
-         ;                     ;   gde3w_n = gde3w_0   !        ---          !
+            gdept_b = gdept_0  ;   gdept_n = gdept_0   !        ---          ! depth of grid-points
+            gdepw_b = gdepw_0  ;   gdepw_n = gdepw_0   !        ---          !
+                                   gde3w_n = gde3w_0   !        ---          !
          !                                                                  
-         ;    e3t_b =   e3t_0  ;     e3t_n =   e3t_0   ;   e3t_a =  e3t_0    ! scale factors
-         ;    e3u_b =   e3u_0  ;     e3u_n =   e3u_0   ;   e3u_a =  e3u_0    !
-         ;    e3v_b =   e3v_0  ;     e3v_n =   e3v_0   ;   e3v_a =  e3v_0    !
-         ;                     ;     e3f_n =   e3f_0   !        ---          !
-         ;    e3w_b =   e3w_0  ;     e3w_n =   e3w_0   !        ---          !
-         ;   e3uw_b =  e3uw_0  ;    e3uw_n =  e3uw_0   !        ---          !
-         ;   e3vw_b =  e3vw_0  ;    e3vw_n =  e3vw_0   !        ---          !
+              e3t_b =   e3t_0  ;     e3t_n =   e3t_0   ;   e3t_a =  e3t_0    ! scale factors
+              e3u_b =   e3u_0  ;     e3u_n =   e3u_0   ;   e3u_a =  e3u_0    !
+              e3v_b =   e3v_0  ;     e3v_n =   e3v_0   ;   e3v_a =  e3v_0    !
+                                     e3f_n =   e3f_0   !        ---          !
+              e3w_b =   e3w_0  ;     e3w_n =   e3w_0   !        ---          !
+             e3uw_b =  e3uw_0  ;    e3uw_n =  e3uw_0   !        ---          !
+             e3vw_b =  e3vw_0  ;    e3vw_n =  e3vw_0   !        ---          !
          !
          z1_hu_0(:,:) = ssumask(:,:) / ( hu_0(:,:) + 1._wp - ssumask(:,:) )     ! _i mask due to ISF
@@ -169,9 +169,9 @@
          !
          !        before       !          now          !       after         !
-         ;                     ;      ht_n =    ht_0   !                     ! water column thickness
-         ;     hu_b =    hu_0  ;      hu_n =    hu_0   ;    hu_a =    hu_0   ! 
-         ;     hv_b =    hv_0  ;      hv_n =    hv_0   ;    hv_a =    hv_0   !
-         ;  r1_hu_b = z1_hu_0  ;   r1_hu_n = z1_hu_0   ; r1_hu_a = z1_hu_0   ! inverse of water column thickness
-         ;  r1_hv_b = z1_hv_0  ;   r1_hv_n = z1_hv_0   ; r1_hv_a = z1_hv_0   !
+                                      ht_n =    ht_0   !                     ! water column thickness
+               hu_b =    hu_0  ;      hu_n =    hu_0   ;    hu_a =    hu_0   ! 
+               hv_b =    hv_0  ;      hv_n =    hv_0   ;    hv_a =    hv_0   !
+            r1_hu_b = z1_hu_0  ;   r1_hu_n = z1_hu_0   ; r1_hu_a = z1_hu_0   ! inverse of water column thickness
+            r1_hv_b = z1_hv_0  ;   r1_hv_n = z1_hv_0   ; r1_hv_a = z1_hv_0   !
          !
          !
Index: branches/2016/dev_CNRS_2016/NEMOGCM/NEMO/OPA_SRC/DYN/dynadv_cen2.F90
===================================================================
--- branches/2016/dev_CNRS_2016/NEMOGCM/NEMO/OPA_SRC/DYN/dynadv_cen2.F90	(revision 7280)
+++ branches/2016/dev_CNRS_2016/NEMOGCM/NEMO/OPA_SRC/DYN/dynadv_cen2.F90	(revision 7282)
@@ -118,6 +118,6 @@
       ENDIF
       DO jk = 2, jpkm1                    ! interior advective fluxes
-         DO jj = 2, jpjm1                       ! 1/4 * Vertical transport
-            DO ji = fs_2, fs_jpim1
+         DO jj = 2, jpj                       ! 1/4 * Vertical transport
+            DO ji = 2, jpi
                zfw(ji,jj,jk) = 0.25_wp * e1e2t(ji,jj) * wn(ji,jj,jk)
             END DO
Index: branches/2016/dev_CNRS_2016/NEMOGCM/NEMO/OPA_SRC/DYN/dynadv_ubs.F90
===================================================================
--- branches/2016/dev_CNRS_2016/NEMOGCM/NEMO/OPA_SRC/DYN/dynadv_ubs.F90	(revision 7280)
+++ branches/2016/dev_CNRS_2016/NEMOGCM/NEMO/OPA_SRC/DYN/dynadv_ubs.F90	(revision 7282)
@@ -211,6 +211,6 @@
       ENDIF
       DO jk = 2, jpkm1                          ! interior fluxes
-         DO jj = 2, jpjm1
-            DO ji = fs_2, fs_jpim1
+         DO jj = 2, jpj
+            DO ji = 2, jpi
                zfw(ji,jj,jk) = 0.25_wp * e1e2t(ji,jj) * wn(ji,jj,jk)
             END DO
Index: branches/2016/dev_CNRS_2016/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg.F90
===================================================================
--- branches/2016/dev_CNRS_2016/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg.F90	(revision 7280)
+++ branches/2016/dev_CNRS_2016/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg.F90	(revision 7282)
@@ -205,5 +205,5 @@
       ENDIF
       !                          ! Control of surface pressure gradient scheme options
-      ;                              nspg =  np_NO    ;   ioptio = 0
+                                     nspg =  np_NO    ;   ioptio = 0
       IF( ln_dynspg_exp ) THEN   ;   nspg =  np_EXP   ;   ioptio = ioptio + 1   ;   ENDIF
       IF( ln_dynspg_ts  ) THEN   ;   nspg =  np_TS    ;   ioptio = ioptio + 1   ;   ENDIF
Index: branches/2016/dev_CNRS_2016/NEMOGCM/NEMO/OPA_SRC/DYN/dynzdf_imp.F90
===================================================================
--- branches/2016/dev_CNRS_2016/NEMOGCM/NEMO/OPA_SRC/DYN/dynzdf_imp.F90	(revision 7280)
+++ branches/2016/dev_CNRS_2016/NEMOGCM/NEMO/OPA_SRC/DYN/dynzdf_imp.F90	(revision 7282)
@@ -294,5 +294,5 @@
             ze3va =  ( 1._wp - r_vvl ) * e3v_n(ji,jj,1) + r_vvl * e3v_a(ji,jj,1) 
             va(ji,jj,1) = va(ji,jj,1) + p2dt * 0.5_wp * ( vtau_b(ji,jj) + vtau(ji,jj) )   &
-               &                                      / ( ze3va * rau0 ) 
+               &                                      / ( ze3va * rau0 ) * vmask(ji,jj,1) 
          END DO
       END DO
Index: branches/2016/dev_CNRS_2016/NEMOGCM/NEMO/OPA_SRC/ICB/icbini.F90
===================================================================
--- branches/2016/dev_CNRS_2016/NEMOGCM/NEMO/OPA_SRC/ICB/icbini.F90	(revision 7280)
+++ branches/2016/dev_CNRS_2016/NEMOGCM/NEMO/OPA_SRC/ICB/icbini.F90	(revision 7282)
@@ -120,5 +120,5 @@
       ! first entry with narea for this processor is left hand interior index
       ! last  entry                               is right hand interior index
-      jj = jpj/2
+      jj = nlcj/2
       nicbdi = -1
       nicbei = -1
@@ -136,5 +136,5 @@
       !
       ! repeat for j direction
-      ji = jpi/2
+      ji = nlci/2
       nicbdj = -1
       nicbej = -1
@@ -153,9 +153,9 @@
       ! special for east-west boundary exchange we save the destination index
       i1 = MAX( nicbdi-1, 1)
-      i3 = INT( src_calving(i1,jpj/2) )
+      i3 = INT( src_calving(i1,nlcj/2) )
       jj = INT( i3/nicbpack )
       ricb_left = REAL( i3 - nicbpack*jj, wp )
       i1 = MIN( nicbei+1, jpi )
-      i3 = INT( src_calving(i1,jpj/2) )
+      i3 = INT( src_calving(i1,nlcj/2) )
       jj = INT( i3/nicbpack )
       ricb_right = REAL( i3 - nicbpack*jj, wp )
@@ -196,5 +196,5 @@
          WRITE(numicb,*) 'berg left       ', ricb_left
          WRITE(numicb,*) 'berg right      ', ricb_right
-         jj = jpj/2
+         jj = nlcj/2
          WRITE(numicb,*) "central j line:"
          WRITE(numicb,*) "i processor"
@@ -202,5 +202,5 @@
          WRITE(numicb,*) "i point"
          WRITE(numicb,*) (INT(src_calving(ji,jj)), ji=1,jpi)
-         ji = jpi/2
+         ji = nlci/2
          WRITE(numicb,*) "central i line:"
          WRITE(numicb,*) "j processor"
Index: branches/2016/dev_CNRS_2016/NEMOGCM/NEMO/OPA_SRC/LBC/lib_mpp.F90
===================================================================
--- branches/2016/dev_CNRS_2016/NEMOGCM/NEMO/OPA_SRC/LBC/lib_mpp.F90	(revision 7280)
+++ branches/2016/dev_CNRS_2016/NEMOGCM/NEMO/OPA_SRC/LBC/lib_mpp.F90	(revision 7282)
@@ -2983,5 +2983,6 @@
       !!----------------------------------------------------------------------
       !
-      ALLOCATE( ztab(jpiglo,4,num_fields), znorthloc(jpi,4,num_fields), zfoldwk(jpi,4,num_fields), znorthgloio(jpi,4,num_fields,jpni) )   ! expanded to 3 dimensions
+      ALLOCATE( ztab(jpiglo,4,num_fields), znorthloc(jpi,4,num_fields), zfoldwk(jpi,4,num_fields),   &
+            &   znorthgloio(jpi,4,num_fields,jpni) )   ! expanded to 3 dimensions
       ALLOCATE( ztabl(jpi,4,num_fields), ztabr(jpi*jpmaxngh, 4,num_fields) )
       !
Index: branches/2016/dev_CNRS_2016/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk.F90
===================================================================
--- branches/2016/dev_CNRS_2016/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk.F90	(revision 7280)
+++ branches/2016/dev_CNRS_2016/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk.F90	(revision 7282)
@@ -281,5 +281,7 @@
       IF( MOD( kt - 1, nn_fsbc ) == 0 )   THEN
          qlw_ice(:,:,1)   = sf(jp_qlw )%fnow(:,:,1)
-         qsr_ice(:,:,1)   = sf(jp_qsr )%fnow(:,:,1)
+         IF( ln_dm2dc ) THEN ; qsr_ice(:,:,1) = sbc_dcy( sf(jp_qsr)%fnow(:,:,1) )
+         ELSE                ; qsr_ice(:,:,1) =          sf(jp_qsr)%fnow(:,:,1) 
+         ENDIF 
          tatm_ice(:,:)    = sf(jp_tair)%fnow(:,:,1)
          qatm_ice(:,:)    = sf(jp_humi)%fnow(:,:,1)
Index: branches/2016/dev_CNRS_2016/NEMOGCM/NEMO/OPA_SRC/SBC/sbccpl.F90
===================================================================
--- branches/2016/dev_CNRS_2016/NEMOGCM/NEMO/OPA_SRC/SBC/sbccpl.F90	(revision 7280)
+++ branches/2016/dev_CNRS_2016/NEMOGCM/NEMO/OPA_SRC/SBC/sbccpl.F90	(revision 7282)
@@ -1327,57 +1327,66 @@
       !!             ***  ROUTINE sbc_cpl_ice_flx  ***
       !!
-      !! ** Purpose :   provide the heat and freshwater fluxes of the 
-      !!              ocean-ice system.
+      !! ** Purpose :   provide the heat and freshwater fluxes of the ocean-ice system
       !!
       !! ** Method  :   transform the fields received from the atmosphere into
       !!             surface heat and fresh water boundary condition for the 
       !!             ice-ocean system. The following fields are provided:
-      !!              * total non solar, solar and freshwater fluxes (qns_tot, 
+      !!               * total non solar, solar and freshwater fluxes (qns_tot, 
       !!             qsr_tot and emp_tot) (total means weighted ice-ocean flux)
       !!             NB: emp_tot include runoffs and calving.
-      !!              * fluxes over ice (qns_ice, qsr_ice, emp_ice) where
+      !!               * fluxes over ice (qns_ice, qsr_ice, emp_ice) where
       !!             emp_ice = sublimation - solid precipitation as liquid
       !!             precipitation are re-routed directly to the ocean and 
-      !!             runoffs and calving directly enter the ocean.
-      !!              * solid precipitation (sprecip), used to add to qns_tot 
+      !!             calving directly enter the ocean (runoffs are read but included in trasbc.F90)
+      !!               * solid precipitation (sprecip), used to add to qns_tot 
       !!             the heat lost associated to melting solid precipitation
       !!             over the ocean fraction.
-      !!       ===>> CAUTION here this changes the net heat flux received from
-      !!             the atmosphere
-      !!
-      !!                  - the fluxes have been separated from the stress as
-      !!                 (a) they are updated at each ice time step compare to
-      !!                 an update at each coupled time step for the stress, and
-      !!                 (b) the conservative computation of the fluxes over the
-      !!                 sea-ice area requires the knowledge of the ice fraction
-      !!                 after the ice advection and before the ice thermodynamics,
-      !!                 so that the stress is updated before the ice dynamics
-      !!                 while the fluxes are updated after it.
+      !!               * heat content of rain, snow and evap can also be provided,
+      !!             otherwise heat flux associated with these mass flux are
+      !!             guessed (qemp_oce, qemp_ice)
+      !!
+      !!             - the fluxes have been separated from the stress as
+      !!               (a) they are updated at each ice time step compare to
+      !!               an update at each coupled time step for the stress, and
+      !!               (b) the conservative computation of the fluxes over the
+      !!               sea-ice area requires the knowledge of the ice fraction
+      !!               after the ice advection and before the ice thermodynamics,
+      !!               so that the stress is updated before the ice dynamics
+      !!               while the fluxes are updated after it.
+      !!
+      !! ** Details
+      !!             qns_tot = pfrld * qns_oce + ( 1 - pfrld ) * qns_ice   => provided
+      !!                     + qemp_oce + qemp_ice                         => recalculated and added up to qns
+      !!
+      !!             qsr_tot = pfrld * qsr_oce + ( 1 - pfrld ) * qsr_ice   => provided
+      !!
+      !!             emp_tot = emp_oce + emp_ice                           => calving is provided and added to emp_tot (and emp_oce)
+      !!                                                                      river runoff (rnf) is provided but not included here
       !!
       !! ** Action  :   update at each nf_ice time step:
       !!                   qns_tot, qsr_tot  non-solar and solar total heat fluxes
       !!                   qns_ice, qsr_ice  non-solar and solar heat fluxes over the ice
-      !!                   emp_tot            total evaporation - precipitation(liquid and solid) (-runoff)(-calving)
-      !!                   emp_ice            ice sublimation - solid precipitation over the ice
-      !!                   dqns_ice           d(non-solar heat flux)/d(Temperature) over the ice
-      !!                   sprecip             solid precipitation over the ocean  
+      !!                   emp_tot           total evaporation - precipitation(liquid and solid) (-calving)
+      !!                   emp_ice           ice sublimation - solid precipitation over the ice
+      !!                   dqns_ice          d(non-solar heat flux)/d(Temperature) over the ice
+      !!                   sprecip           solid precipitation over the ocean  
       !!----------------------------------------------------------------------
-      REAL(wp), INTENT(in   ), DIMENSION(:,:)             ::   p_frld  ! lead fraction            [0 to 1]
+      REAL(wp), INTENT(in   ), DIMENSION(:,:)   ::   p_frld     ! lead fraction                [0 to 1]
       ! optional arguments, used only in 'mixed oce-ice' case
-      REAL(wp), INTENT(in   ), DIMENSION(:,:,:), OPTIONAL ::   palbi   ! all skies ice albedo 
-      REAL(wp), INTENT(in   ), DIMENSION(:,:  ), OPTIONAL ::   psst    ! sea surface temperature  [Celsius]
-      REAL(wp), INTENT(in   ), DIMENSION(:,:,:), OPTIONAL ::   pist    ! ice surface temperature  [Kelvin]
-      !
-      INTEGER ::   jl   ! dummy loop index
+      REAL(wp), INTENT(in   ), DIMENSION(:,:,:), OPTIONAL ::   palbi      ! all skies ice albedo 
+      REAL(wp), INTENT(in   ), DIMENSION(:,:  ), OPTIONAL ::   psst       ! sea surface temperature     [Celsius]
+      REAL(wp), INTENT(in   ), DIMENSION(:,:,:), OPTIONAL ::   pist       ! ice surface temperature     [Kelvin]
+      !
+      INTEGER ::   jl         ! dummy loop index
       REAL(wp), POINTER, DIMENSION(:,:  ) ::   zcptn, ztmp, zicefr, zmsk, zsnw
-      REAL(wp), POINTER, DIMENSION(:,:  ) ::   zemp_tot, zemp_ice, zemp_oce, ztprecip, zsprecip, zevap, zevap_ice, zdevap_ice
+      REAL(wp), POINTER, DIMENSION(:,:  ) ::   zemp_tot, zemp_ice, zemp_oce, ztprecip, zsprecip, zevap_oce, zevap_ice, zdevap_ice
       REAL(wp), POINTER, DIMENSION(:,:  ) ::   zqns_tot, zqns_oce, zqsr_tot, zqsr_oce, zqprec_ice, zqemp_oce, zqemp_ice
       REAL(wp), POINTER, DIMENSION(:,:,:) ::   zqns_ice, zqsr_ice, zdqns_ice, zqevap_ice
       !!----------------------------------------------------------------------
       !
-      IF( nn_timing == 1 )   CALL timing_start('sbc_cpl_ice_flx')
+      IF( nn_timing == 1 )  CALL timing_start('sbc_cpl_ice_flx')
       !
       CALL wrk_alloc( jpi,jpj,     zcptn, ztmp, zicefr, zmsk, zsnw )
-      CALL wrk_alloc( jpi,jpj,     zemp_tot, zemp_ice, zemp_oce, ztprecip, zsprecip, zevap, zevap_ice, zdevap_ice )
+      CALL wrk_alloc( jpi,jpj,     zemp_tot, zemp_ice, zemp_oce, ztprecip, zsprecip, zevap_oce, zevap_ice, zdevap_ice )
       CALL wrk_alloc( jpi,jpj,     zqns_tot, zqns_oce, zqsr_tot, zqsr_oce, zqprec_ice, zqemp_oce, zqemp_ice )
       CALL wrk_alloc( jpi,jpj,jpl, zqns_ice, zqsr_ice, zdqns_ice, zqevap_ice )
@@ -1388,39 +1397,41 @@
       !
       !                                                      ! ========================= !
-      !                                                      !    freshwater budget      !   (emp)
+      !                                                      !    freshwater budget      !   (emp_tot)
       !                                                      ! ========================= !
       !
-      !                                                           ! total Precipitation - total Evaporation (emp_tot)
-      !                                                           ! solid precipitation - sublimation       (emp_ice)
-      !                                                           ! solid Precipitation                     (sprecip)
-      !                                                           ! liquid + solid Precipitation            (tprecip)
+      !                                                           ! solid Precipitation                                (sprecip)
+      !                                                           ! liquid + solid Precipitation                       (tprecip)
+      !                                                           ! total Evaporation - total Precipitation            (emp_tot)
+      !                                                           ! sublimation - solid precipitation (cell average)   (emp_ice)
       SELECT CASE( TRIM( sn_rcv_emp%cldes ) )
-      CASE( 'conservative'  )   ! received fields: jpr_rain, jpr_snow, jpr_ievp, jpr_tevp
-         zsprecip(:,:) = frcv(jpr_snow)%z3(:,:,1)                  ! May need to ensure positive here
-         ztprecip(:,:) = frcv(jpr_rain)%z3(:,:,1) + zsprecip(:,:)  ! May need to ensure positive here
-         zemp_tot(:,:) = frcv(jpr_tevp)%z3(:,:,1) - ztprecip(:,:)
-         zemp_ice(:,:) = frcv(jpr_ievp)%z3(:,:,1) - frcv(jpr_snow)%z3(:,:,1)
-            CALL iom_put( 'rain'         , frcv(jpr_rain)%z3(:,:,1)              )   ! liquid precipitation 
+      CASE( 'conservative' )   ! received fields: jpr_rain, jpr_snow, jpr_ievp, jpr_tevp
+         zsprecip(:,:) =   frcv(jpr_snow)%z3(:,:,1)                  ! May need to ensure positive here
+         ztprecip(:,:) =   frcv(jpr_rain)%z3(:,:,1) + zsprecip(:,:)  ! May need to ensure positive here
+         zemp_tot(:,:) =   frcv(jpr_tevp)%z3(:,:,1) - ztprecip(:,:)
+         zemp_ice(:,:) = ( frcv(jpr_ievp)%z3(:,:,1) - frcv(jpr_snow)%z3(:,:,1) ) * zicefr(:,:)
+               CALL iom_put( 'rain'         ,   frcv(jpr_rain)%z3(:,:,1)                                                         )  ! liquid precipitation 
          IF( iom_use('hflx_rain_cea') )   &
-            CALL iom_put( 'hflx_rain_cea', frcv(jpr_rain)%z3(:,:,1) * zcptn(:,:) )   ! heat flux from liq. precip. 
-         IF( iom_use('evap_ao_cea') .OR. iom_use('hflx_evap_cea') )   &
-            ztmp(:,:) = frcv(jpr_tevp)%z3(:,:,1) - frcv(jpr_ievp)%z3(:,:,1) * zicefr(:,:)
+            &  CALL iom_put( 'hflx_rain_cea',   frcv(jpr_rain)%z3(:,:,1) * zcptn(:,:)                                            )  ! heat flux from liq. precip. 
          IF( iom_use('evap_ao_cea'  ) )   &
-            CALL iom_put( 'evap_ao_cea'  , ztmp                   )   ! ice-free oce evap (cell average)
+            &  CALL iom_put( 'evap_ao_cea'  ,   frcv(jpr_tevp)%z3(:,:,1) - frcv(jpr_ievp)%z3(:,:,1) * zicefr(:,:)                )  ! ice-free oce evap (cell average)
          IF( iom_use('hflx_evap_cea') )   &
-            CALL iom_put( 'hflx_evap_cea', ztmp(:,:) * zcptn(:,:) )   ! heat flux from from evap (cell average)
-      CASE( 'oce and ice'   )   ! received fields: jpr_sbpr, jpr_semp, jpr_oemp, jpr_ievp
+            &  CALL iom_put( 'hflx_evap_cea', ( frcv(jpr_tevp)%z3(:,:,1) - frcv(jpr_ievp)%z3(:,:,1) * zicefr(:,:) ) * zcptn(:,:) )  ! heat flux from from evap (cell average)
+      CASE( 'oce and ice' )   ! received fields: jpr_sbpr, jpr_semp, jpr_oemp, jpr_ievp
          zemp_tot(:,:) = p_frld(:,:) * frcv(jpr_oemp)%z3(:,:,1) + zicefr(:,:) * frcv(jpr_sbpr)%z3(:,:,1)
-         zemp_ice(:,:) = frcv(jpr_semp)%z3(:,:,1)
+         zemp_ice(:,:) = frcv(jpr_semp)%z3(:,:,1) * zicefr(:,:)
          zsprecip(:,:) = frcv(jpr_ievp)%z3(:,:,1) - frcv(jpr_semp)%z3(:,:,1)
          ztprecip(:,:) = frcv(jpr_semp)%z3(:,:,1) - frcv(jpr_sbpr)%z3(:,:,1) + zsprecip(:,:)
       END SELECT
+
 #if defined key_lim3
-      ! zsnw = snow percentage over ice after wind blowing
-      zsnw(:,:) = 0._wp
-      CALL lim_thd_snwblow( p_frld, zsnw )
+      ! zsnw = snow fraction over ice after wind blowing
+      zsnw(:,:) = 0._wp  ;  CALL lim_thd_snwblow( p_frld, zsnw )
       
-      ! --- evaporation (used later in sbccpl) --- !
-      zevap(:,:) = zemp_tot(:,:) + ztprecip(:,:)
+      ! --- evaporation minus precipitation corrected (because of wind blowing on snow) --- !
+      zemp_ice(:,:) = zemp_ice(:,:) + zsprecip(:,:) * ( zicefr(:,:) - zsnw(:,:) )  ! emp_ice = A * sublimation - zsnw * sprecip
+      zemp_oce(:,:) = zemp_tot(:,:) - zemp_ice(:,:)                                ! emp_oce = emp_tot - emp_ice
+
+      ! --- evaporation over ocean (used later for qemp) --- !
+      zevap_oce(:,:) = frcv(jpr_tevp)%z3(:,:,1) - frcv(jpr_ievp)%z3(:,:,1) * zicefr(:,:)
 
       ! --- evaporation over ice (kg/m2/s) --- !
@@ -1430,8 +1441,4 @@
       zdevap_ice(:,:) = 0._wp
       
-      ! --- evaporation minus precipitation corrected for the effect of wind blowing on snow --- !
-      zemp_oce(:,:) = zemp_tot(:,:) - zemp_ice(:,:) - zsprecip * (1._wp - zsnw)
-      zemp_ice(:,:) = zemp_ice(:,:) + zsprecip * (1._wp - zsnw)    
-
       ! --- runoffs (included in emp later on) --- !
       IF( srcv(jpr_rnf)%laction )   rnf(:,:) = frcv(jpr_rnf)%z3(:,:,1)
@@ -1443,5 +1450,5 @@
          CALL iom_put( 'calving_cea', frcv(jpr_cal)%z3(:,:,1) )
       ENDIF
-      
+
       IF( ln_mixcpl ) THEN
          emp_tot(:,:) = emp_tot(:,:) * xcplmask(:,:,0) + zemp_tot(:,:) * zmsk(:,:)
@@ -1465,12 +1472,10 @@
          ENDDO
       ENDIF
-      
-      IF( iom_use('subl_ai_cea') )   CALL iom_put( 'subl_ai_cea', zevap_ice(:,:) * zicefr(:,:)    )  ! Sublimation over sea-ice (cell average)
-                                     CALL iom_put( 'snowpre'    , sprecip                         )  ! Snow
-      IF( iom_use('snow_ao_cea') )   CALL iom_put( 'snow_ao_cea', sprecip(:,:) * ( 1._wp - zsnw ) )  ! Snow over ice-free ocean  (cell average)
-      IF( iom_use('snow_ai_cea') )   CALL iom_put( 'snow_ai_cea', sprecip(:,:) *           zsnw   )  ! Snow over sea-ice         (cell average)    
+
+      IF( iom_use('subl_ai_cea') )   CALL iom_put( 'subl_ai_cea', zevap_ice(:,:) * zicefr(:,:)         )  ! Sublimation over sea-ice (cell average)
+                                     CALL iom_put( 'snowpre'    , sprecip(:,:)                         )  ! Snow
+      IF( iom_use('snow_ao_cea') )   CALL iom_put( 'snow_ao_cea', sprecip(:,:) * ( 1._wp - zsnw(:,:) ) )  ! Snow over ice-free ocean  (cell average)
+      IF( iom_use('snow_ai_cea') )   CALL iom_put( 'snow_ai_cea', sprecip(:,:) *           zsnw(:,:)   )  ! Snow over sea-ice         (cell average)
 #else
-      ! Sublimation over sea-ice (cell average)
-      IF( iom_use('subl_ai_cea') )  CALL iom_put( 'subl_ai_cea', frcv(jpr_ievp)%z3(:,:,1) * zicefr(:,:) )
       ! runoffs and calving (put in emp_tot)
       IF( srcv(jpr_rnf)%laction )   rnf(:,:) = frcv(jpr_rnf)%z3(:,:,1)
@@ -1492,9 +1497,8 @@
       ENDIF
 
-         CALL iom_put( 'snowpre'    , sprecip                                )   ! Snow
-      IF( iom_use('snow_ao_cea') )   &
-         CALL iom_put( 'snow_ao_cea', sprecip(:,:) * p_frld(:,:)             )   ! Snow        over ice-free ocean  (cell average)
-      IF( iom_use('snow_ai_cea') )   &
-         CALL iom_put( 'snow_ai_cea', sprecip(:,:) * zicefr(:,:)             )   ! Snow        over sea-ice         (cell average)
+      IF( iom_use('subl_ai_cea') )  CALL iom_put( 'subl_ai_cea', frcv(jpr_ievp)%z3(:,:,1) * zicefr(:,:) )  ! Sublimation over sea-ice (cell average)
+                                    CALL iom_put( 'snowpre'    , sprecip(:,:)               )   ! Snow
+      IF( iom_use('snow_ao_cea') )  CALL iom_put( 'snow_ao_cea', sprecip(:,:) * p_frld(:,:) )   ! Snow over ice-free ocean  (cell average)
+      IF( iom_use('snow_ai_cea') )  CALL iom_put( 'snow_ai_cea', sprecip(:,:) * zicefr(:,:) )   ! Snow over sea-ice         (cell average)
 #endif
 
@@ -1502,18 +1506,17 @@
       SELECT CASE( TRIM( sn_rcv_qns%cldes ) )                !   non solar heat fluxes   !   (qns)
       !                                                      ! ========================= !
-      CASE( 'oce only' )                                     ! the required field is directly provided
-         zqns_tot(:,:  ) = frcv(jpr_qnsoce)%z3(:,:,1)
-      CASE( 'conservative' )                                      ! the required fields are directly provided
-         zqns_tot(:,:  ) = frcv(jpr_qnsmix)%z3(:,:,1)
+      CASE( 'oce only' )         ! the required field is directly provided
+         zqns_tot(:,:) = frcv(jpr_qnsoce)%z3(:,:,1)
+      CASE( 'conservative' )     ! the required fields are directly provided
+         zqns_tot(:,:) = frcv(jpr_qnsmix)%z3(:,:,1)
          IF ( TRIM(sn_rcv_qns%clcat) == 'yes' ) THEN
             zqns_ice(:,:,1:jpl) = frcv(jpr_qnsice)%z3(:,:,1:jpl)
          ELSE
-            ! Set all category values equal for the moment
             DO jl=1,jpl
-               zqns_ice(:,:,jl) = frcv(jpr_qnsice)%z3(:,:,1)
+               zqns_ice(:,:,jl) = frcv(jpr_qnsice)%z3(:,:,1) ! Set all category values equal
             ENDDO
          ENDIF
-      CASE( 'oce and ice' )       ! the total flux is computed from ocean and ice fluxes
-         zqns_tot(:,:  ) =  p_frld(:,:) * frcv(jpr_qnsoce)%z3(:,:,1)
+      CASE( 'oce and ice' )      ! the total flux is computed from ocean and ice fluxes
+         zqns_tot(:,:) =  p_frld(:,:) * frcv(jpr_qnsoce)%z3(:,:,1)
          IF ( TRIM(sn_rcv_qns%clcat) == 'yes' ) THEN
             DO jl=1,jpl
@@ -1522,5 +1525,5 @@
             ENDDO
          ELSE
-            qns_tot(:,:   ) = qns_tot(:,:) + zicefr(:,:) * frcv(jpr_qnsice)%z3(:,:,1)
+            qns_tot(:,:) = qns_tot(:,:) + zicefr(:,:) * frcv(jpr_qnsice)%z3(:,:,1)
             DO jl=1,jpl
                zqns_tot(:,:   ) = zqns_tot(:,:) + zicefr(:,:) * frcv(jpr_qnsice)%z3(:,:,1)
@@ -1528,10 +1531,10 @@
             ENDDO
          ENDIF
-      CASE( 'mixed oce-ice' )     ! the ice flux is cumputed from the total flux, the SST and ice informations
+      CASE( 'mixed oce-ice' )    ! the ice flux is cumputed from the total flux, the SST and ice informations
 ! ** NEED TO SORT OUT HOW THIS SHOULD WORK IN THE MULTI-CATEGORY CASE - CURRENTLY NOT ALLOWED WHEN INTERFACE INITIALISED **
          zqns_tot(:,:  ) = frcv(jpr_qnsmix)%z3(:,:,1)
          zqns_ice(:,:,1) = frcv(jpr_qnsmix)%z3(:,:,1)    &
             &            + frcv(jpr_dqnsdt)%z3(:,:,1) * ( pist(:,:,1) - ( (rt0 + psst(:,:  ) ) * p_frld(:,:)   &
-            &                                                   +          pist(:,:,1)   * zicefr(:,:) ) )
+            &                                           + pist(:,:,1) * zicefr(:,:) ) )
       END SELECT
 !!gm
@@ -1543,13 +1546,9 @@
 !! similar job should be done for snow and precipitation temperature
       !                                     
-      IF( srcv(jpr_cal)%laction ) THEN                            ! Iceberg melting 
-         ztmp(:,:) = frcv(jpr_cal)%z3(:,:,1) * lfus               ! add the latent heat of iceberg melting 
-         zqns_tot(:,:) = zqns_tot(:,:) - ztmp(:,:)
-         IF( iom_use('hflx_cal_cea') )   &
-            CALL iom_put( 'hflx_cal_cea', ztmp + frcv(jpr_cal)%z3(:,:,1) * zcptn(:,:) )   ! heat flux from calving
-      ENDIF
-
-      ztmp(:,:) = p_frld(:,:) * zsprecip(:,:) * lfus
-      IF( iom_use('hflx_snow_cea') )    CALL iom_put( 'hflx_snow_cea', ztmp + sprecip(:,:) * zcptn(:,:) )   ! heat flux from snow (cell average)
+      IF( srcv(jpr_cal)%laction ) THEN   ! Iceberg melting 
+         zqns_tot(:,:) = zqns_tot(:,:) - frcv(jpr_cal)%z3(:,:,1) * lfus  ! add the latent heat of iceberg melting
+                                                                         ! we suppose it melts at 0deg, though it should be temp. of surrounding ocean
+         IF( iom_use('hflx_cal_cea') )   CALL iom_put( 'hflx_cal_cea', - frcv(jpr_cal)%z3(:,:,1) * lfus )   ! heat flux from calving
+      ENDIF
 
 #if defined key_lim3      
@@ -1560,18 +1559,18 @@
 
       ! --- heat flux associated with emp (W/m2) --- !
-      zqemp_oce(:,:) = -      zevap(:,:)                   * p_frld(:,:)      *   zcptn(:,:)   &      ! evap
-         &             + ( ztprecip(:,:) - zsprecip(:,:) )                    *   zcptn(:,:)   &      ! liquid precip
-         &             +   zsprecip(:,:)                   * ( 1._wp - zsnw ) * ( zcptn(:,:) - lfus ) ! solid precip over ocean
+      zqemp_oce(:,:) = -  zevap_oce(:,:)                                      *   zcptn(:,:)   &       ! evap
+         &             + ( ztprecip(:,:) - zsprecip(:,:) )                    *   zcptn(:,:)   &       ! liquid precip
+         &             +   zsprecip(:,:)                   * ( 1._wp - zsnw ) * ( zcptn(:,:) - lfus )  ! solid precip over ocean + snow melting
 !      zqemp_ice(:,:) = -   frcv(jpr_ievp)%z3(:,:,1)        * zicefr(:,:)      *   zcptn(:,:)   &      ! ice evap
 !         &             +   zsprecip(:,:)                   * zsnw             * ( zcptn(:,:) - lfus ) ! solid precip over ice
       zqemp_ice(:,:) =      zsprecip(:,:)                   * zsnw             * ( zcptn(:,:) - lfus ) ! solid precip over ice (only)
-                                                                                                       ! qevap_ice=0 since we consider Tice=0°C
+                                                                                                       ! qevap_ice=0 since we consider Tice=0degC
       
-      ! --- heat content of precip over ice in J/m3 (to be used in 1D-thermo) --- !
+      ! --- enthalpy of snow precip over ice in J/m3 (to be used in 1D-thermo) --- !
       zqprec_ice(:,:) = rhosn * ( zcptn(:,:) - lfus )
 
       ! --- heat content of evap over ice in W/m2 (to be used in 1D-thermo) --- !
       DO jl = 1, jpl
-         zqevap_ice(:,:,jl) = 0._wp ! should be -evap * ( ( Tice - rt0 ) * cpic ) but we do not have Tice, so we consider Tice=0°C
+         zqevap_ice(:,:,jl) = 0._wp ! should be -evap * ( ( Tice - rt0 ) * cpic ) but we do not have Tice, so we consider Tice=0degC
       END DO
 
@@ -1599,10 +1598,18 @@
          qemp_ice (:,:  ) = zqemp_ice (:,:  )
       ENDIF
+
+      !! clem: we should output qemp_oce and qemp_ice (at least)
+      IF( iom_use('hflx_snow_cea') )   CALL iom_put( 'hflx_snow_cea', sprecip(:,:) * ( zcptn(:,:) - Lfus ) )   ! heat flux from snow (cell average)
+      !! these diags are not outputed yet
+!!      IF( iom_use('hflx_rain_cea') )   CALL iom_put( 'hflx_rain_cea', ( tprecip(:,:) - sprecip(:,:) ) * zcptn(:,:) )   ! heat flux from rain (cell average)
+!!      IF( iom_use('hflx_snow_ao_cea') ) CALL iom_put( 'hflx_snow_ao_cea', sprecip(:,:) * ( zcptn(:,:) - Lfus ) * (1._wp - zsnw(:,:)) ) ! heat flux from snow (cell average)
+!!      IF( iom_use('hflx_snow_ai_cea') ) CALL iom_put( 'hflx_snow_ai_cea', sprecip(:,:) * ( zcptn(:,:) - Lfus ) * zsnw(:,:) ) ! heat flux from snow (cell average)
+
 #else
-      ! clem: this formulation is certainly wrong... but better than it was before...
+      ! clem: this formulation is certainly wrong... but better than it was...
       zqns_tot(:,:) = zqns_tot(:,:)                       &            ! zqns_tot update over free ocean with:
          &          - ztmp(:,:)                           &            ! remove the latent heat flux of solid precip. melting
          &          - (  zemp_tot(:,:)                    &            ! remove the heat content of mass flux (assumed to be at SST)
-         &             - zemp_ice(:,:) * zicefr(:,:)  ) * zcptn(:,:) 
+         &             - zemp_ice(:,:) ) * zcptn(:,:) 
 
      IF( ln_mixcpl ) THEN
@@ -1616,6 +1623,6 @@
          qns_ice(:,:,:) = zqns_ice(:,:,:)
       ENDIF
-      !
 #endif
+
       !                                                      ! ========================= !
       SELECT CASE( TRIM( sn_rcv_qsr%cldes ) )                !      solar heat fluxes    !   (qsr)
@@ -1724,5 +1731,5 @@
 
       CALL wrk_dealloc( jpi,jpj,     zcptn, ztmp, zicefr, zmsk, zsnw )
-      CALL wrk_dealloc( jpi,jpj,     zemp_tot, zemp_ice, zemp_oce, ztprecip, zsprecip, zevap, zevap_ice, zdevap_ice )
+      CALL wrk_dealloc( jpi,jpj,     zemp_tot, zemp_ice, zemp_oce, ztprecip, zsprecip, zevap_oce, zevap_ice, zdevap_ice )
       CALL wrk_dealloc( jpi,jpj,     zqns_tot, zqns_oce, zqsr_tot, zqsr_oce, zqprec_ice, zqemp_oce, zqemp_ice )
       CALL wrk_dealloc( jpi,jpj,jpl, zqns_ice, zqsr_ice, zdqns_ice, zqevap_ice )
Index: branches/2016/dev_CNRS_2016/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfric.F90
===================================================================
--- branches/2016/dev_CNRS_2016/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfric.F90	(revision 7280)
+++ branches/2016/dev_CNRS_2016/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfric.F90	(revision 7282)
@@ -162,5 +162,8 @@
                   &                               + avmv(ji,jj,jk) + avmv(ji,jj-1,jk)  )   &
                   &          + avtb(jk) * tmask(ji,jj,jk)
-               !                                            ! Add the background coefficient on eddy viscosity
+            END DO
+         END DO
+         DO jj = 2, jpjm1                                   ! Add the background coefficient on eddy viscosity
+            DO ji = fs_2, fs_jpim1
                avmu(ji,jj,jk) = avmu(ji,jj,jk) + avmb(jk) * umask(ji,jj,jk)
                avmv(ji,jj,jk) = avmv(ji,jj,jk) + avmb(jk) * vmask(ji,jj,jk)
Index: branches/2016/dev_CNRS_2016/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zche.F90
===================================================================
--- branches/2016/dev_CNRS_2016/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zche.F90	(revision 7280)
+++ branches/2016/dev_CNRS_2016/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zche.F90	(revision 7282)
@@ -31,6 +31,7 @@
    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(:,:,:) ::   chemc    ! Solubilities of O2 and CO2
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   chemo2   ! Solubilities of O2 and CO2
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   tempis   ! In situ temperature
 
    REAL(wp), PUBLIC ::   atcox  = 0.20946         ! units atm
@@ -39,38 +40,9 @@
    REAL(wp) ::   o2atm  = 1. / ( 1000. * 0.20946 )  
 
-   REAL(wp) ::   akcc1  = -171.9065       ! coeff. for apparent solubility equilibrium
-   REAL(wp) ::   akcc2  =   -0.077993     ! Millero et al. 1995 from Mucci 1983
-   REAL(wp) ::   akcc3  = 2839.319        
-   REAL(wp) ::   akcc4  =   71.595        
-   REAL(wp) ::   akcc5  =   -0.77712      
-   REAL(wp) ::   akcc6  =    0.00284263   
-   REAL(wp) ::   akcc7  =  178.34        
-   REAL(wp) ::   akcc8  =   -0.07711     
-   REAL(wp) ::   akcc9  =    0.0041249   
-
-   REAL(wp) ::   rgas   = 83.143         ! universal gas constants
+   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) ::   ca0    = -162.8301      ! WEISS & PRICE 1980, units mol/(kg atm)
-   REAL(wp) ::   ca1    =  218.2968
-   REAL(wp) ::   ca2    =   90.9241
-   REAL(wp) ::   ca3    =   -1.47696
-   REAL(wp) ::   ca4    =    0.025695
-   REAL(wp) ::   ca5    =   -0.025225
-   REAL(wp) ::   ca6    =    0.0049867
-
-   REAL(wp) ::   c10    = -3670.7        ! Coeff. for 1. dissoc. of carbonic acid (Edmond and Gieskes, 1970)   
-   REAL(wp) ::   c11    =    62.008     
-   REAL(wp) ::   c12    =    -9.7944    
-   REAL(wp) ::   c13    =     0.0118     
-   REAL(wp) ::   c14    =    -0.000116
-
-   REAL(wp) ::   c20    = -1394.7       ! coeff. for 2. dissoc. of carbonic acid (Millero, 1995)   
-   REAL(wp) ::   c21    =    -4.777   
-   REAL(wp) ::   c22    =     0.0184   
-   REAL(wp) ::   c23    =    -0.000118
 
    REAL(wp) ::   st1    =      0.14     ! constants for calculate concentrations for sulfate
@@ -144,6 +116,6 @@
       REAL(wp) ::   ztgg , ztgg2, ztgg3 , ztgg4 , ztgg5
       REAL(wp) ::   zpres, ztc  , zcl   , zcpexp, zoxy  , zcpexp2
-      REAL(wp) ::   zsqrt, ztr  , zlogt , zcek1
-      REAL(wp) ::   zis  , zis2 , zsal15, zisqrt
+      REAL(wp) ::   zsqrt, ztr  , zlogt , zcek1, zc1, zplat
+      REAL(wp) ::   zis  , zis2 , zsal15, zisqrt, za1  , za2
       REAL(wp) ::   zckb , zck1 , zck2  , zckw  , zak1 , zak2  , zakb , zaksp0, zakw
       REAL(wp) ::   zst  , zft  , zcks  , zckf  , zaksp1
@@ -151,4 +123,20 @@
       !
       IF( nn_timing == 1 )  CALL timing_start('p4z_che')
+      !
+      ! Computations of chemical constants require in situ temperature
+      ! Here a quite simple formulation is used to convert 
+      ! potential temperature to in situ temperature. The errors is less than 
+      ! 0.04°C relative to an exact computation
+      ! ---------------------------------------------------------------------
+      DO jk = 1, jpk
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               zpres = gdept_n(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) )
+               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
+            END DO
+         END DO
+      END DO
       !
       ! CHEMICAL CONSTANTS - SURFACE LAYER
@@ -157,5 +145,5 @@
          DO ji = 1, jpi
             !                             ! SET ABSOLUTE TEMPERATURE
-            ztkel = tsn(ji,jj,1,jp_tem) + 273.15
+            ztkel = tempis(ji,jj,1) + 273.15
             zt    = ztkel * 0.01
             zt2   = zt * zt
@@ -165,7 +153,10 @@
             !                             ! LN(K0) OF SOLUBILITY OF CO2 (EQ. 12, WEISS, 1980)
             !                             !     AND FOR THE ATMOSPHERE FOR NON IDEAL GAS
-            zcek1 = ca0 + ca1 / zt + ca2 * zlogt + ca3 * zt2 + zsal * ( ca4 + ca5 * zt + ca6 * zt2 )
+            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) = EXP( zcek1 ) * 1.e-6 * rhop(ji,jj,1) / 1000.  ! mol/(L uatm)
+            chemc(ji,jj,1) = EXP( zcek1 ) * 1.e-6 * rhop(ji,jj,1) / 1000. ! mol/(kg uatm)
+            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
             !
          END DO
@@ -177,8 +168,8 @@
          DO jj = 1, jpj
             DO ji = 1, jpi
-              ztkel = tsn(ji,jj,jk,jp_tem) + 273.15
+              ztkel = tempis(ji,jj,jk) + 273.15
               zsal  = tsn(ji,jj,jk,jp_sal) + ( 1.- tmask(ji,jj,jk) ) * 35.
               zsal2 = zsal * zsal
-              ztgg  = LOG( ( 298.15 - tsn(ji,jj,jk,jp_tem) ) / ztkel )  ! Set the GORDON & GARCIA scaled temperature
+              ztgg  = LOG( ( 298.15 - tempis(ji,jj,jk) ) / ztkel )  ! Set the GORDON & GARCIA scaled temperature
               ztgg2 = ztgg  * ztgg
               ztgg3 = ztgg2 * ztgg
@@ -200,9 +191,12 @@
             DO ji = 1, jpi
 
-               ! SET PRESSION
-               zpres   = 1.025e-1 * gdept_n(ji,jj,jk)
+               ! SET PRESSION ACCORDING TO SAUNDER (1980)
+               zplat   = SIN ( ABS(gphit(ji,jj)*3.141592654/180.) )
+               zc1 = 5.92E-3 + zplat**2 * 5.25E-3
+               zpres = ((1-zc1)-SQRT(((1-zc1)**2)-(8.84E-6*gdept_n(ji,jj,jk)))) / 4.42E-6
+               zpres = zpres / 10.0
 
                ! SET ABSOLUTE TEMPERATURE
-               ztkel   = tsn(ji,jj,jk,jp_tem) + 273.15
+               ztkel   = tempis(ji,jj,jk) + 273.15
                zsal    = tsn(ji,jj,jk,jp_sal) + ( 1.-tmask(ji,jj,jk) ) * 35.
                zsqrt  = SQRT( zsal )
@@ -213,5 +207,5 @@
                zis2   = zis * zis
                zisqrt = SQRT( zis )
-               ztc     = tsn(ji,jj,jk,jp_tem) + ( 1.- tmask(ji,jj,jk) ) * 20.
+               ztc     = tempis(ji,jj,jk) + ( 1.- tmask(ji,jj,jk) ) * 20.
 
                ! CHLORINITY (WOOSTER ET AL., 1969)
@@ -246,6 +240,10 @@
 
 
-               zck1    = c10 * ztr + c11 + c12 * zlogt + c13 * zsal + c14 * zsal * zsal
-               zck2    = c20 * ztr + c21 + c22 * zsal   + c23 * zsal**2
+               ! DISSOCIATION COEFFICIENT FOR CARBONATE ACCORDING TO 
+               ! MEHRBACH (1973) REFIT BY MILLERO (1995), seawater scale
+               zck1    = -1.0*(3633.86*ztr - 61.2172 + 9.6777*zlogt  &
+                  - 0.011555*zsal + 0.0001152*zsal*zsal)
+               zck2    = -1.0*(471.78*ztr + 25.9290 - 3.16967*zlogt      &
+                  - 0.01781*zsal + 0.0001122*zsal*zsal)
 
                ! PKW (H2O) (DICKSON AND RILEY, 1979)
@@ -256,6 +254,7 @@
                ! APPARENT SOLUBILITY PRODUCT K'SP OF CALCITE IN SEAWATER
                !       (S=27-43, T=2-25 DEG C) at pres =0 (atmos. pressure) (MUCCI 1983)
-               zaksp0  = akcc1 + akcc2 * ztkel + akcc3 * ztr + akcc4 * LOG10( ztkel )   &
-                  &   + ( akcc5 + akcc6 * ztkel + akcc7 * ztr ) * zsqrt + akcc8 * zsal + akcc9 * zsal15
+               zaksp0  = -171.9065 -0.077993*ztkel + 2839.319*ztr + 71.595*LOG10( ztkel )   &
+                  &      + (-0.77712 + 0.00284263*ztkel + 178.34*ztr) * zsqrt  &
+                  &      - 0.07711*zsal + 0.0041249*zsal15
 
                ! K1, K2 OF CARBONIC ACID, KB OF BORIC ACID, KW (H2O) (LIT.?)
@@ -327,6 +326,6 @@
       !!                     ***  ROUTINE p4z_che_alloc  ***
       !!----------------------------------------------------------------------
-      ALLOCATE( sio3eq(jpi,jpj,jpk), fekeq(jpi,jpj,jpk), chemc(jpi,jpj), chemo2(jpi,jpj,jpk),   &
-      &         STAT=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 )
       !
       IF( p4z_che_alloc /= 0 )   CALL ctl_warn('p4z_che_alloc : failed to allocate arrays.')
Index: branches/2016/dev_CNRS_2016/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zflx.F90
===================================================================
--- branches/2016/dev_CNRS_2016/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zflx.F90	(revision 7280)
+++ branches/2016/dev_CNRS_2016/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zflx.F90	(revision 7282)
@@ -84,13 +84,14 @@
       REAL(wp) ::   ztc, ztc2, ztc3, ztc4, zws, zkgwan
       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) ::   zyr_dec, zdco2dt
       CHARACTER (len=25) :: charout
-      REAL(wp), POINTER, DIMENSION(:,:) :: zkgco2, zkgo2, zh2co3, zoflx, zw2d 
+      REAL(wp), POINTER, DIMENSION(:,:) :: zkgco2, zkgo2, zh2co3, zoflx, zw2d, zpco2atm  
       !!---------------------------------------------------------------------
       !
       IF( nn_timing == 1 )  CALL timing_start('p4z_flx')
       !
-      CALL wrk_alloc( jpi, jpj, zkgco2, zkgo2, zh2co3, zoflx )
+      CALL wrk_alloc( jpi, jpj, zkgco2, zkgo2, zh2co3, zoflx, zpco2atm )
       !
 
@@ -177,15 +178,24 @@
       DO jj = 1, jpj
          DO ji = 1, jpi
+            ztkel     = tsn(ji,jj,1,jp_tem) + 273.15
+            zsal      = tsn(ji,jj,1,jp_sal) + ( 1.- tmask(ji,jj,1) ) * 35.
+            zvapsw    = EXP(24.4543 - 67.4509*(100.0/ztkel) - 4.8489*LOG(ztkel/100) - 0.000544*zsal)
+            zpco2atm(ji,jj) = satmco2(ji,jj) * ( patm(ji,jj) - zvapsw )
+            zxc2      = ( 1.0 - zpco2atm(ji,jj) * 1E-6 )**2
+            zfugcoeff = EXP( patm(ji,jj) * (chemc(ji,jj,2) + 2.0 * zxc2 * chemc(ji,jj,3) )   &
+            &           / ( 82.05736 * ztkel ))
+            zfco2 = zpco2atm(ji,jj) * zfugcoeff
+
             ! Compute CO2 flux for the sea and air
-            zfld = satmco2(ji,jj) * patm(ji,jj) * tmask(ji,jj,1) * chemc(ji,jj) * zkgco2(ji,jj)   ! (mol/L) * (m/s)
-            zflu = zh2co3(ji,jj) * tmask(ji,jj,1) * zkgco2(ji,jj)                                   ! (mol/L) (m/s) ?
+            zfld = zfco2 * chemc(ji,jj,1) * zkgco2(ji,jj)  ! (mol/L) * (m/s)
+            zflu = zh2co3(ji,jj) * zkgco2(ji,jj)                                   ! (mol/L) (m/s) ?
             oce_co2(ji,jj) = ( zfld - zflu ) * rfact2 * e1e2t(ji,jj) * tmask(ji,jj,1) * 1000.
             ! compute the trend
-            tra(ji,jj,1,jpdic) = tra(ji,jj,1,jpdic) + ( zfld - zflu ) * rfact2 / e3t_n(ji,jj,1)
+            tra(ji,jj,1,jpdic) = tra(ji,jj,1,jpdic) + ( zfld - zflu ) * rfact2 / e3t_n(ji,jj,1) * tmask(ji,jj,1)
 
             ! Compute O2 flux 
-            zfld16 = patm(ji,jj) * chemo2(ji,jj,1) * tmask(ji,jj,1) * zkgo2(ji,jj)          ! (mol/L) * (m/s)
-            zflu16 = trb(ji,jj,1,jpoxy) * tmask(ji,jj,1) * zkgo2(ji,jj)
-            zoflx(ji,jj) = zfld16 - zflu16
+            zfld16 = patm(ji,jj) * chemo2(ji,jj,1) * zkgo2(ji,jj)          ! (mol/L) * (m/s)
+            zflu16 = trb(ji,jj,1,jpoxy) * zkgo2(ji,jj)
+            zoflx(ji,jj) = ( zfld16 - zflu16 ) * tmask(ji,jj,1)
             tra(ji,jj,1,jpoxy) = tra(ji,jj,1,jpoxy) + zoflx(ji,jj) * rfact2 / e3t_n(ji,jj,1)
          END DO
@@ -218,5 +228,5 @@
          ENDIF
          IF( iom_use( "Dpco2" ) ) THEN
-           zw2d(:,:) = ( satmco2(:,:) * patm(:,:) - zh2co3(:,:) / ( chemc(:,:) + rtrn ) ) * tmask(:,:,1)
+           zw2d(:,:) = ( zpco2atm(:,:) - zh2co3(:,:) / ( chemc(:,:,1) + rtrn ) ) * tmask(:,:,1)
            CALL iom_put( "Dpco2" ,  zw2d )
          ENDIF
@@ -234,9 +244,9 @@
             trc2d(:,:,jp_pcs0_2d + 1) = zoflx(:,:) * 1000 * tmask(:,:,1) 
             trc2d(:,:,jp_pcs0_2d + 2) = zkgco2(:,:) * tmask(:,:,1) 
-            trc2d(:,:,jp_pcs0_2d + 3) = ( satmco2(:,:) * patm(:,:) - zh2co3(:,:) / ( chemc(:,:) + rtrn ) ) * tmask(:,:,1) 
-         ENDIF
-      ENDIF
-      !
-      CALL wrk_dealloc( jpi, jpj, zkgco2, zkgo2, zh2co3, zoflx )
+            trc2d(:,:,jp_pcs0_2d + 3) = ( zpco2atm(:,:) - zh2co3(:,:) / ( chemc(:,:,1) + rtrn ) ) * tmask(:,:,1) 
+         ENDIF
+      ENDIF
+      !
+      CALL wrk_dealloc( jpi, jpj, zkgco2, zkgo2, zh2co3, zoflx, zpco2atm )
       !
       IF( nn_timing == 1 )  CALL timing_stop('p4z_flx')
Index: branches/2016/dev_CNRS_2016/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zlim.F90
===================================================================
--- branches/2016/dev_CNRS_2016/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zlim.F90	(revision 7280)
+++ branches/2016/dev_CNRS_2016/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zlim.F90	(revision 7282)
@@ -44,4 +44,5 @@
    REAL(wp), PUBLIC ::  xkdoc       !:  2nd half-sat. of DOC remineralization  
    REAL(wp), PUBLIC ::  concbfe     !:  Fe half saturation for bacteria 
+   REAL(wp), PUBLIC ::  oxymin      !:  half saturation constant for anoxia
    REAL(wp), PUBLIC ::  qnfelim     !:  optimal Fe quota for nanophyto
    REAL(wp), PUBLIC ::  qdfelim     !:  optimal Fe quota for diatoms
@@ -186,4 +187,14 @@
       END DO
       !
+      DO jk = 1, jpkm1
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               ! denitrification factor computed from O2 levels
+               nitrfac(ji,jj,jk) = MAX(  0.e0, 0.4 * ( 6.e-6  - trb(ji,jj,jk,jpoxy) )    &
+                  &                                / ( oxymin + trb(ji,jj,jk,jpoxy) )  )
+               nitrfac(ji,jj,jk) = MIN( 1., nitrfac(ji,jj,jk) )
+            END DO
+         END DO
+      END DO
       !
       IF( lk_iomput .AND. knt == nrdttrc ) THEN        ! save output diagnostics
@@ -215,5 +226,5 @@
       NAMELIST/nampislim/ concnno3, concdno3, concnnh4, concdnh4, concnfer, concdfer, concbfe,   &
          &                concbno3, concbnh4, xsizedia, xsizephy, xsizern, xsizerd,          & 
-         &                xksi1, xksi2, xkdoc, qnfelim, qdfelim, caco3r
+         &                xksi1, xksi2, xkdoc, qnfelim, qdfelim, caco3r, oxymin
       INTEGER :: ios                 ! Local integer output status for namelist read
 
@@ -248,8 +259,11 @@
          WRITE(numout,*) '    Minimum size criteria for nanophyto      xsizephy  = ', xsizephy
          WRITE(numout,*) '    Fe half saturation for bacteria          concbfe   = ', concbfe
+         WRITE(numout,*) '    halk saturation constant for anoxia       oxymin   =' , oxymin
          WRITE(numout,*) '    optimal Fe quota for nano.               qnfelim   = ', qnfelim
          WRITE(numout,*) '    Optimal Fe quota for diatoms             qdfelim   = ', qdfelim
       ENDIF
-
+      !
+      nitrfac (:,:,:) = 0._wp
+      !
    END SUBROUTINE p4z_lim_init
 
Index: branches/2016/dev_CNRS_2016/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zlys.F90
===================================================================
--- branches/2016/dev_CNRS_2016/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zlys.F90	(revision 7280)
+++ branches/2016/dev_CNRS_2016/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zlys.F90	(revision 7282)
@@ -65,10 +65,10 @@
       REAL(wp) ::   zomegaca, zexcess, zexcess0
       CHARACTER (len=25) :: charout
-      REAL(wp), POINTER, DIMENSION(:,:,:) :: zco3, zcaldiss   
+      REAL(wp), POINTER, DIMENSION(:,:,:) :: zco3, zco3sat, zcaldiss   
       !!---------------------------------------------------------------------
       !
       IF( nn_timing == 1 )  CALL timing_start('p4z_lys')
       !
-      CALL wrk_alloc( jpi, jpj, jpk, zco3, zcaldiss )
+      CALL wrk_alloc( jpi, jpj, jpk, zco3, zco3sat, zcaldiss )
       !
       zco3    (:,:,:) = 0.
@@ -117,5 +117,6 @@
                zcalcon  = calcon * ( tsn(ji,jj,jk,jp_sal) / 35._wp )
                zfact    = rhop(ji,jj,jk) / 1000._wp
-               zomegaca = ( zcalcon * zco3(ji,jj,jk) * zfact ) / aksp(ji,jj,jk) 
+               zomegaca = ( zcalcon * zco3(ji,jj,jk) ) / ( aksp(ji,jj,jk) * zfact + rtrn )
+               zco3sat(ji,jj,jk) = aksp(ji,jj,jk) * zfact / ( zcalcon + rtrn )
 
                ! SET DEGREE OF UNDER-/SUPERSATURATION
@@ -146,12 +147,12 @@
       IF( lk_iomput .AND. knt == nrdttrc ) THEN
          IF( iom_use( "PH"     ) ) CALL iom_put( "PH"    , -1. * LOG10( hi(:,:,:) )          * tmask(:,:,:) )
-         IF( iom_use( "CO3"    ) ) CALL iom_put( "CO3"   , zco3(:,:,:) * 1.e+3               * tmask(:,:,:) )
-         IF( iom_use( "CO3sat" ) ) CALL iom_put( "CO3sat", aksp(:,:,:) * 1.e+3 / calcon      * tmask(:,:,:) )
-         IF( iom_use( "DCAL"   ) ) CALL iom_put( "DCAL"  , zcaldiss(:,:,:) * 1.e+3 * rfact2r   * tmask(:,:,:) )
+         IF( iom_use( "CO3"    ) ) CALL iom_put( "CO3"   , zco3(:,:,:)    * 1.e+3            * tmask(:,:,:) )
+         IF( iom_use( "CO3sat" ) ) CALL iom_put( "CO3sat", zco3sat(:,:,:) * 1.e+3            * tmask(:,:,:) )
+         IF( iom_use( "DCAL"   ) ) CALL iom_put( "DCAL"  , zcaldiss(:,:,:) * 1.e+3 * rfact2r * tmask(:,:,:) )
       ELSE
          IF( ln_diatrc ) THEN
             trc3d(:,:,:,jp_pcs0_3d    ) = -1. * LOG10( hi(:,:,:) ) * tmask(:,:,:)
             trc3d(:,:,:,jp_pcs0_3d + 1) = zco3(:,:,:)              * tmask(:,:,:)
-            trc3d(:,:,:,jp_pcs0_3d + 2) = aksp(:,:,:) / calcon     * tmask(:,:,:)
+            trc3d(:,:,:,jp_pcs0_3d + 2) = zco3sat(:,:,:)           * tmask(:,:,:)
          ENDIF
       ENDIF
@@ -163,5 +164,5 @@
       ENDIF
       !
-      CALL wrk_dealloc( jpi, jpj, jpk, zco3, zcaldiss )
+      CALL wrk_dealloc( jpi, jpj, jpk, zco3, zco3sat, zcaldiss )
       !
       IF( nn_timing == 1 )  CALL timing_stop('p4z_lys')
Index: branches/2016/dev_CNRS_2016/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zopt.F90
===================================================================
--- branches/2016/dev_CNRS_2016/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zopt.F90	(revision 7280)
+++ branches/2016/dev_CNRS_2016/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zopt.F90	(revision 7282)
@@ -110,10 +110,8 @@
       !                                        !  --------------------------------------
       IF( l_trcdm2dc ) THEN                     !  diurnal cycle
-         !                                       ! 1% of qsr to compute euphotic layer
-         zqsr100(:,:)   = 0.01 * qsr_mean(:,:)     !  daily mean qsr
          !
          zqsr_corr(:,:) = qsr_mean(:,:) / ( 1. - fr_i(:,:) + rtrn )
          !
-         CALL p4z_opt_par( kt, zqsr_corr, ze1, ze2, ze3 ) 
+         CALL p4z_opt_par( kt, zqsr_corr, ze1, ze2, ze3, pqsr100 = zqsr100 ) 
          !
          DO jk = 1, nksrp      
@@ -132,10 +130,8 @@
          !
       ELSE
-         !                                       ! 1% of qsr to compute euphotic layer
-         zqsr100(:,:)   = 0.01 * qsr(:,:)     !  daily mean qsr
          !
          zqsr_corr(:,:) = qsr(:,:) / ( 1. - fr_i(:,:) + rtrn )
          !
-         CALL p4z_opt_par( kt, zqsr_corr, ze1, ze2, ze3 ) 
+         CALL p4z_opt_par( kt, zqsr_corr, ze1, ze2, ze3, pqsr100 = zqsr100  ) 
          !
          DO jk = 1, nksrp      
@@ -165,5 +161,5 @@
          DO jj = 1, jpj
            DO ji = 1, jpi
-              IF( etot_ndcy(ji,jj,jk) * tmask(ji,jj,jk) >= 0.43 * zqsr100(ji,jj) )  THEN
+              IF( etot_ndcy(ji,jj,jk) * tmask(ji,jj,jk) >=  zqsr100(ji,jj) )  THEN
                  neln(ji,jj) = jk+1                    ! Euphotic level : 1rst T-level strictly below Euphotic layer
                  !                                     ! nb: ensure the compatibility with nmld_trc definition in trd_mld_trc_zint
@@ -234,5 +230,5 @@
    END SUBROUTINE p4z_opt
 
-   SUBROUTINE p4z_opt_par( kt, pqsr, pe1, pe2, pe3, pe0 ) 
+   SUBROUTINE p4z_opt_par( kt, pqsr, pe1, pe2, pe3, pe0, pqsr100 ) 
       !!----------------------------------------------------------------------
       !!                  ***  routine p4z_opt_par  ***
@@ -247,4 +243,5 @@
       REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout)           ::  pe1 , pe2 , pe3   !  PAR ( R-G-B)
       REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout), OPTIONAL ::  pe0  
+      REAL(wp), DIMENSION(jpi,jpj)    , INTENT(out)  , OPTIONAL  ::  pqsr100  
       !! * local variables
       INTEGER    ::   ji, jj, jk     ! dummy loop indices
@@ -256,5 +253,8 @@
       ELSE                  ;  zqsr(:,:) = xparsw         * pqsr(:,:)
       ENDIF
-      !
+      
+      !  Light at the euphotic depth 
+      IF( PRESENT( pqsr100 ) )  pqsr100(:,:) = 0.01 * 3. * zqsr(:,:)
+
       IF( PRESENT( pe0 ) ) THEN     !  W-level
          !
Index: branches/2016/dev_CNRS_2016/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zprod.F90
===================================================================
--- branches/2016/dev_CNRS_2016/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zprod.F90	(revision 7280)
+++ branches/2016/dev_CNRS_2016/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zprod.F90	(revision 7282)
@@ -279,4 +279,6 @@
                   zprdia(ji,jj,jk) = zprdia(ji,jj,jk) * zmixdiat(ji,jj)
                ENDIF
+               zprbio(ji,jj,jk) = zprbio(ji,jj,jk) * ( 1. - fr_i(ji,jj) )
+               zprdia(ji,jj,jk) = zprdia(ji,jj,jk) * ( 1. - fr_i(ji,jj) )
             END DO
          END DO
@@ -315,51 +317,28 @@
       END DO
 
-      IF( ln_newprod ) THEN
-         DO jk = 1, jpkm1
-            DO jj = 1, jpj
-               DO ji = 1, jpi
-                  IF( gdepw_n(ji,jj,jk+1) <= hmld(ji,jj) ) THEN
-                     zprnch(ji,jj,jk) = zprnch(ji,jj,jk) * zmixnano(ji,jj)
-                     zprdch(ji,jj,jk) = zprdch(ji,jj,jk) * zmixdiat(ji,jj)
-                  ENDIF
-                  IF( etot_ndcy(ji,jj,jk) > 1.E-3 ) THEN
-                     !  production terms for nanophyto. ( chlorophyll )
-                     znanotot = enano(ji,jj,jk) * zstrn(ji,jj)
-                     zprod    = rday * zprorca(ji,jj,jk) * zprnch(ji,jj,jk) * xlimphy(ji,jj,jk)
-                     zprochln(ji,jj,jk) = chlcmin * 12. * zprorca (ji,jj,jk)
-                     zprochln(ji,jj,jk) = zprochln(ji,jj,jk) + (chlcnm-chlcmin) * 12. * zprod / &
-                                        & (  zpislopead(ji,jj,jk) * znanotot +rtrn)
-                     !  production terms for diatomees ( chlorophyll )
-                     zdiattot = ediat(ji,jj,jk) * zstrn(ji,jj)
-                     zprod = rday * zprorcad(ji,jj,jk) * zprdch(ji,jj,jk) * xlimdia(ji,jj,jk)
-                     zprochld(ji,jj,jk) = chlcmin * 12. * zprorcad(ji,jj,jk)
-                     zprochld(ji,jj,jk) = zprochld(ji,jj,jk) + (chlcdm-chlcmin) * 12. * zprod / &
-                                        & ( zpislopead2(ji,jj,jk) * zdiattot +rtrn )
-                  ENDIF
-               END DO
-            END DO
-         END DO
-      ELSE
-         DO jk = 1, jpkm1
-            DO jj = 1, jpj
-               DO ji = 1, jpi
-                  IF( etot_ndcy(ji,jj,jk) > 1.E-3 ) THEN
-                     !  production terms for nanophyto. ( chlorophyll )
-                     znanotot = enano(ji,jj,jk)
-                     zprod = rday * zprorca(ji,jj,jk) * zprnch(ji,jj,jk) * trb(ji,jj,jk,jpphy) * xlimphy(ji,jj,jk)
-                     zprochln(ji,jj,jk) = chlcmin * 12. * zprorca (ji,jj,jk)
-                     zprochln(ji,jj,jk) = zprochln(ji,jj,jk) + (chlcnm-chlcmin) * 144. * zprod            &
-                     &                    / ( zpislopead(ji,jj,jk) * trb(ji,jj,jk,jpnch) * znanotot +rtrn )
-                     !  production terms for diatomees ( chlorophyll )
-                     zdiattot = ediat(ji,jj,jk)
-                     zprod = rday * zprorcad(ji,jj,jk) * zprdch(ji,jj,jk) * trb(ji,jj,jk,jpdia) * xlimdia(ji,jj,jk)
-                     zprochld(ji,jj,jk) = chlcmin * 12. * zprorcad(ji,jj,jk)
-                     zprochld(ji,jj,jk) = zprochld(ji,jj,jk) + (chlcdm-chlcmin) * 144. * zprod             &
-                     &                    / ( zpislopead2(ji,jj,jk) * trb(ji,jj,jk,jpdch) * zdiattot +rtrn )
-                  ENDIF
-               END DO
-            END DO
-         END DO
-      ENDIF
+      DO jk = 1, jpkm1
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               IF( gdepw_n(ji,jj,jk+1) <= hmld(ji,jj) ) THEN
+                  zprnch(ji,jj,jk) = zprnch(ji,jj,jk) * zmixnano(ji,jj)
+                  zprdch(ji,jj,jk) = zprdch(ji,jj,jk) * zmixdiat(ji,jj)
+               ENDIF
+               IF( etot_ndcy(ji,jj,jk) > 1.E-3 ) THEN
+                  !  production terms for nanophyto. ( chlorophyll )
+                  znanotot = enano(ji,jj,jk) * zstrn(ji,jj)
+                  zprod    = rday * zprorca(ji,jj,jk) * zprnch(ji,jj,jk) * xlimphy(ji,jj,jk)
+                  zprochln(ji,jj,jk) = chlcmin * 12. * zprorca (ji,jj,jk)
+                  zprochln(ji,jj,jk) = zprochln(ji,jj,jk) + (chlcnm-chlcmin) * 12. * zprod / &
+                                     & (  zpislopead(ji,jj,jk) * znanotot +rtrn)
+                  !  production terms for diatomees ( chlorophyll )
+                  zdiattot = ediat(ji,jj,jk) * zstrn(ji,jj)
+                  zprod = rday * zprorcad(ji,jj,jk) * zprdch(ji,jj,jk) * xlimdia(ji,jj,jk)
+                  zprochld(ji,jj,jk) = chlcmin * 12. * zprorcad(ji,jj,jk)
+                  zprochld(ji,jj,jk) = zprochld(ji,jj,jk) + (chlcdm-chlcmin) * 12. * zprod / &
+                                     & ( zpislopead2(ji,jj,jk) * zdiattot +rtrn )
+               ENDIF
+            END DO
+         END DO
+      END DO
 
       !   Update the arrays TRA which contain the biological sources and sinks
Index: branches/2016/dev_CNRS_2016/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zrem.F90
===================================================================
--- branches/2016/dev_CNRS_2016/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zrem.F90	(revision 7280)
+++ branches/2016/dev_CNRS_2016/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zrem.F90	(revision 7282)
@@ -44,5 +44,4 @@
    REAL(wp), PUBLIC ::  xsiremlab  !: fast remineralisation rate of POC 
    REAL(wp), PUBLIC ::  xsilab     !: fraction of labile biogenic silica 
-   REAL(wp), PUBLIC ::  oxymin     !: halk saturation constant for anoxia 
 
 
@@ -109,15 +108,4 @@
                   zdepprod(ji,jj,jk) = zdepmin**0.273
                ENDIF
-            END DO
-         END DO
-      END DO
-
-      DO jk = 1, jpkm1
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-               ! denitrification factor computed from O2 levels
-               nitrfac(ji,jj,jk) = MAX(  0.e0, 0.4 * ( 6.e-6  - trb(ji,jj,jk,jpoxy) )    &
-                  &                                / ( oxymin + trb(ji,jj,jk,jpoxy) )  )
-               nitrfac(ji,jj,jk) = MIN( 1., nitrfac(ji,jj,jk) )
             END DO
          END DO
@@ -355,6 +343,5 @@
       !!
       !!----------------------------------------------------------------------
-      NAMELIST/nampisrem/ xremik, xremip, nitrif, xsirem, xsiremlab, xsilab,   &
-      &                   oxymin
+      NAMELIST/nampisrem/ xremik, xremip, nitrif, xsirem, xsiremlab, xsilab
       INTEGER :: ios                 ! Local integer output status for namelist read
 
@@ -378,8 +365,6 @@
          WRITE(numout,*) '    fraction of labile biogenic silica        xsilab    =', xsilab
          WRITE(numout,*) '    NH4 nitrification rate                    nitrif    =', nitrif
-         WRITE(numout,*) '    halk saturation constant for anoxia       oxymin    =', oxymin
       ENDIF
       !
-      nitrfac (:,:,:) = 0._wp
       denitr  (:,:,:) = 0._wp
       denitnh4(:,:,:) = 0._wp
Index: branches/2016/dev_CNRS_2016/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zsbc.F90
===================================================================
--- branches/2016/dev_CNRS_2016/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zsbc.F90	(revision 7280)
+++ branches/2016/dev_CNRS_2016/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zsbc.F90	(revision 7282)
@@ -155,10 +155,11 @@
       IF( ln_ndepo ) THEN
          IF( kt == nit000 .OR. ( kt /= nit000 .AND. ntimes_ndep > 1 ) ) THEN
-            CALL fld_read( kt, 1, sf_ndepo )
-            DO jj = 1, jpj
-               DO ji = 1, jpi
-                  nitdep(ji,jj) = sf_ndepo(1)%fnow(ji,jj,1) / rno3 / ( 14E6 * ryyss * e3t_n(ji,jj,1) + rtrn )
-               END DO
-            END DO
+             zcoef = rno3 * 14E6 * ryyss
+             CALL fld_read( kt, 1, sf_ndepo )
+             nitdep(:,:) = sf_ndepo(1)%fnow(:,:,1) / zcoef / e3t_n(:,:,1) 
+         ENDIF
+         IF( .NOT.ln_linssh ) THEN
+           zcoef = rno3 * 14E6 * ryyss
+           nitdep(:,:) = sf_ndepo(1)%fnow(:,:,1) / zcoef / e3t_n(:,:,1) 
          ENDIF
       ENDIF
@@ -461,5 +462,5 @@
          ironsed(:,:,jpk) = 0._wp
          DO jk = 1, jpkm1
-            ironsed(:,:,jk) = sedfeinput * zcmask(:,:,jk) / ( e3t_n(:,:,jk) * rday )
+            ironsed(:,:,jk) = sedfeinput * zcmask(:,:,jk) / ( e3t_0(:,:,jk) * rday )
          END DO
          DEALLOCATE( zcmask)
@@ -479,5 +480,7 @@
          CALL iom_close( numhydro )
          !
-         hydrofe(:,:,:) = ( hydrofe(:,:,:) * hratio ) / ( cvol(:,:,:) * ryyss + rtrn ) / 1000._wp
+         DO jk = 1, jpk
+            hydrofe(:,:,jk) = ( hydrofe(:,:,jk) * hratio ) / ( e1e2t(:,:) * e3t_0(:,:,jk) * ryyss + rtrn ) / 1000._wp
+         ENDDO
          !
       ENDIF
Index: branches/2016/dev_CNRS_2016/NEMOGCM/NEMO/TOP_SRC/TRP/trcsbc.F90
===================================================================
--- branches/2016/dev_CNRS_2016/NEMOGCM/NEMO/TOP_SRC/TRP/trcsbc.F90	(revision 7280)
+++ branches/2016/dev_CNRS_2016/NEMOGCM/NEMO/TOP_SRC/TRP/trcsbc.F90	(revision 7282)
@@ -100,5 +100,5 @@
          IF(lwp) WRITE(numout,*) '~~~~~~~ '
 
-         IF( ln_rsttr .AND. .NOT.ln_top_euler .AND.    &                     ! Restart: read in restart  file
+         IF( ln_rsttr .AND. .NOT.ln_top_euler .AND.   &                     ! Restart: read in restart  file
             iom_varid( numrtr, 'sbc_'//TRIM(ctrcnm(1))//'_b', ldstop = .FALSE. ) > 0 ) THEN
             IF(lwp) WRITE(numout,*) '          nittrc000-nn_dttrc surface tracer content forcing fields red in the restart file'
Index: branches/2016/dev_CNRS_2016/NEMOGCM/NEMO/TOP_SRC/trcrst.F90
===================================================================
--- branches/2016/dev_CNRS_2016/NEMOGCM/NEMO/TOP_SRC/trcrst.F90	(revision 7280)
+++ branches/2016/dev_CNRS_2016/NEMOGCM/NEMO/TOP_SRC/trcrst.F90	(revision 7282)
@@ -295,5 +295,5 @@
          IF(lwp) WRITE(numout,9000) jn, TRIM( ctrcnm(jn) ), zmean, zmin, zmax, zdrift
       END DO
-      WRITE(numout,*) 
+      IF(lwp) WRITE(numout,*) 
 9000  FORMAT(' tracer nb :',i2,'    name :',a10,'    mean :',e18.10,'    min :',e18.10, &
       &      '    max :',e18.10,'    drift :',e18.10, ' %')
Index: branches/2016/dev_CNRS_2016/NEMOGCM/NEMO/TOP_SRC/trcstp.F90
===================================================================
--- branches/2016/dev_CNRS_2016/NEMOGCM/NEMO/TOP_SRC/trcstp.F90	(revision 7280)
+++ branches/2016/dev_CNRS_2016/NEMOGCM/NEMO/TOP_SRC/trcstp.F90	(revision 7282)
@@ -32,6 +32,6 @@
    REAL(wp), DIMENSION(:,:,:), SAVE, ALLOCATABLE ::   qsr_arr ! save qsr during TOP time-step
    REAL(wp) :: rdt_sampl
-   INTEGER  :: nb_rec_per_days
-   INTEGER  :: isecfst, iseclast
+   INTEGER  :: nb_rec_per_day
+   REAL(wp) :: rsecfst, rseclast
    LOGICAL  :: llnew
 
@@ -108,10 +108,9 @@
       END DO
       IF( lwp ) WRITE(numstr,9300) kt,  ztrai / areatot
-9300  FORMAT(i10,e18.10)
+9300  FORMAT(i10,D23.16)
       !
       IF( nn_timing == 1 )   CALL timing_stop('trc_stp')
       !
    END SUBROUTINE trc_stp
-
 
    SUBROUTINE trc_mean_qsr( kt )
@@ -122,5 +121,5 @@
       !!               of diurnal cycle
       !!
-      !! ** Method  : store in TOP the qsr every hour ( or every time-step the latter 
+      !! ** Method  : store in TOP the qsr every hour ( or every time-step if the latter 
       !!              is greater than 1 hour ) and then, compute the  mean with 
       !!              a moving average over 24 hours. 
@@ -129,43 +128,84 @@
       INTEGER, INTENT(in) ::   kt
       INTEGER  :: jn
-      !!----------------------------------------------------------------------
-      !
+      REAL(wp) :: zkt
+      CHARACTER(len=1)               ::   cl1                      ! 1 character
+      CHARACTER(len=2)               ::   cl2                      ! 2 characters
+
       IF( kt == nittrc000 ) THEN
          IF( ln_cpl )  THEN  
-            rdt_sampl = 86400. / ncpl_qsr_freq
-            nb_rec_per_days = ncpl_qsr_freq
+            rdt_sampl = rday / ncpl_qsr_freq
+            nb_rec_per_day = ncpl_qsr_freq
          ELSE  
-            rdt_sampl = MAX( 3600., rdt * nn_dttrc )
-            nb_rec_per_days = INT( 86400 / rdt_sampl )
+            rdt_sampl = MAX( 3600., rdttrc )
+            nb_rec_per_day = INT( rday / rdt_sampl )
          ENDIF
          !
          IF( lwp ) THEN
             WRITE(numout,*) 
-            WRITE(numout,*) ' Sampling frequency dt = ', rdt_sampl, 's','   Number of sampling per day  nrec = ', nb_rec_per_days
+            WRITE(numout,*) ' Sampling frequency dt = ', rdt_sampl, 's','   Number of sampling per day  nrec = ', nb_rec_per_day
             WRITE(numout,*) 
          ENDIF
          !
-         ALLOCATE( qsr_arr(jpi,jpj,nb_rec_per_days ) )
-         DO jn = 1, nb_rec_per_days
-            qsr_arr(:,:,jn) = qsr(:,:)
+         ALLOCATE( qsr_arr(jpi,jpj,nb_rec_per_day ) )
+         !
+         !                                            !* 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
+            rsecfst = INT( zkt ) * rdttrc
+            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
+         ELSE                                         !* no restart: set from nit000 values
+            IF(lwp) WRITE(numout,*) 'trc_qsr_mean:   qsr_mean set to nit000 values'
+            rsecfst  = kt * rdttrc
+            !
+            qsr_mean(:,:) = qsr(:,:)
+            DO jn = 1, nb_rec_per_day
+               qsr_arr(:,:,jn) = qsr_mean(:,:)
+            ENDDO
+         ENDIF
+         !
+      ENDIF
+      !
+      rseclast = kt * rdttrc
+      !
+      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, &
+             &                      ' time = ', rseclast/3600.,'hours '
+          rsecfst = rseclast
+          DO jn = 1, nb_rec_per_day - 1
+             qsr_arr(:,:,jn) = qsr_arr(:,:,jn+1)
+          ENDDO
+          qsr_arr (:,:,nb_rec_per_day) = qsr(:,:)
+          qsr_mean(:,:                ) = SUM( qsr_arr(:,:,:), 3 ) / nb_rec_per_day
+      ENDIF
+      !
+      IF( lrst_trc ) THEN    !* Write the mean of qsr in restart file 
+         IF(lwp) WRITE(numout,*)
+         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 )
+          DO jn = 1, nb_rec_per_day 
+             IF( jn <= 9 )  THEN
+               WRITE(cl1,'(i1)') jn
+               CALL iom_rstput( kt, nitrst, numrtw, 'qsr_arr_'//cl1, qsr_arr(:,:,jn) )
+             ELSE
+               WRITE(cl2,'(i2.2)') jn
+               CALL iom_rstput( kt, nitrst, numrtw, 'qsr_arr_'//cl2, qsr_arr(:,:,jn) )
+             ENDIF
          ENDDO
-         qsr_mean(:,:) = qsr(:,:)
-         !
-         isecfst  = nsec_year + nsec1jan000   !   number of seconds between Jan. 1st 00h of nit000 year and the middle of time step
-         iseclast = isecfst
-         !
-      ENDIF
-      !
-      iseclast = nsec_year + nsec1jan000
-      llnew   = ( iseclast - isecfst )  > INT( rdt_sampl )   !   new shortwave to store
-      IF( kt /= nittrc000 .AND. llnew ) THEN
-          IF( lwp ) WRITE(numout,*) ' New shortwave to sample for TOP at time kt = ', kt, &
-             &                      ' time = ', (iseclast+rdt*nn_dttrc/2.)/3600.,'hours '
-          isecfst = iseclast
-          DO jn = 1, nb_rec_per_days - 1
-             qsr_arr(:,:,jn) = qsr_arr(:,:,jn+1)
-          END DO
-          qsr_arr (:,:,nb_rec_per_days) = qsr(:,:)
-          qsr_mean(:,:                ) = SUM( qsr_arr(:,:,:), 3 ) / nb_rec_per_days
+         CALL iom_rstput( kt, nitrst, numrtw, 'qsr_mean', qsr_mean(:,:) )
       ENDIF
       !
Index: branches/2016/dev_CNRS_2016/NEMOGCM/NEMO/TOP_SRC/trcsub.F90
===================================================================
--- branches/2016/dev_CNRS_2016/NEMOGCM/NEMO/TOP_SRC/trcsub.F90	(revision 7280)
+++ branches/2016/dev_CNRS_2016/NEMOGCM/NEMO/TOP_SRC/trcsub.F90	(revision 7282)
@@ -16,10 +16,4 @@
    USE in_out_manager
    USE lbclnk
-#if defined key_zdftke
-   USE zdftke          ! twice TKE (en)
-#endif
-#if defined key_zdfgls
-   USE zdfgls, ONLY: en
-#endif
    USE trabbl
    USE zdf_oce
Index: branches/2016/dev_CNRS_2016/NEMOGCM/TOOLS/MISCELLANEOUS/README_uspcfg
===================================================================
--- branches/2016/dev_CNRS_2016/NEMOGCM/TOOLS/MISCELLANEOUS/README_uspcfg	(revision 7282)
+++ branches/2016/dev_CNRS_2016/NEMOGCM/TOOLS/MISCELLANEOUS/README_uspcfg	(revision 7282)
@@ -0,0 +1,134 @@
+Instructions for using unsupported configurations on systems without wget access to
+the internet
+
+The unsupported configurations in NEMO are provided as a way for users to provide
+alternative configurations to the wider community without having to rely on system team
+sponsership and support. The idea is to provide minimal but sufficient information
+with the standard distribution to allow ancillary files to be fetched from remote
+servers at compile time. Thus a makenemo command such as:
+
+./makenemo -n MYISOMIP -u ISOMIP -m some_arch
+
+will create a new configuration called MYISOMIP with MY_SRC and EXP00 directories
+populated from a remote server. It does this by first using the wget command to fetch a
+list of remote files and then using wget repeatedly to fetch each of the files listed
+therein. The location of the remote list is contained in one line of information held
+in the uspcfg.txt file (in CONFIG).
+
+These instructions provide a work-around for systems without direct access to the
+internet or on which wget has been disabled/not installed. On these systems it is
+possible to create local copies of the remote files and replace the wget function
+with a local copy command. It will be the user's responsibility to ensure that the
+local copies reflect any changes made in the remote sources.
+
+The procedure requires the following steps:
+
+1. From the target NEMOGCM directory on your target system copy the following
+   files into a working directory on a system that does have internet access
+   and a working wget command:
+
+      a. TOOLS/MISCELLANEOUS/make_usp_tar.sh
+      b. CONFIG/uspcfg.txt
+
+2. On that alternative system and in the working directory containing the newly
+   copied files, run the make_usp_tar.sh script. The script takes three arguments:
+
+      a. The name of the uspcfg.txt file
+      b. The name of the unsupported configuration you wish to use
+      c. The name of a directory that will be created and filled with the remote
+         files in the current working directory. This directory will be tarred and
+         eventually transferred to the target system. There it will be unpacked and
+         will form the local archive for this configuration. A meaningful name is
+         therefore recommended.
+
+   For example:
+
+   ./make_usp_tar.sh uspcfg.txt ISOMIP ISOMIP_ARCHIVE
+
+3. Copy the resulting tarball (in this case: ISOMIP_ARCHIVE.tar) to the target 
+   system and unpack (tar xvf ISOMIP_ARCHIVE.tar) in a location that is visible to 
+   the node on which you intend to compile nemo. The unpacked directory will contain
+   the files retrieved from the remote server, modified versions of the file list and
+   uspcfg.txt files and two new scripts. For example:
+
+   ls -1R ISOMIP_ARCHIVE
+   cpp_ISOMIP.fcm             - downloaded from remote server
+   def_wget                   - newly constructed script
+   local_template.ctl         - file list modified for local copy
+   uspcfg_local_template.txt  - a copy of uspcfg.txt with modification for local copy
+   remote_file.list           - downloaded from remote server
+   set_local_uspcfg           - newly constructed script
+   EXP00                      - new subdirectory
+   MY_SRC                     - new subdirectory
+
+   ./EXP00:
+   iodef.xml                  - downloaded from remote server
+   namelist_cfg               - downloaded from remote server
+
+   ./MY_SRC:
+   domzgr.F90                 - downloaded from remote server
+   istate.F90                 - downloaded from remote server
+
+4. The final stage is to enter the directory and execute the set_local_uspcfg script.
+   This script will prompt for a full path to the target CONFIG directory and then:
+
+   a. Edit the templates and create local.ctl and uspcfg_local.txt files by inserting 
+      the full path to the archive directory (i.e. the current working directory)
+   b. move uspcfg_local.txt to the named CONFIG directory
+   c. rename any existing uspcfg.txt file in the CONFIG directory to uspcfg_remote.txt
+   d. insert symbolic link (uspcfg.txt) to uspcfg_local.txt in the CONFIG directory
+   e. redefine/define wget as an exported bash function to simply copy
+
+   This last action is needed so that subsequent use of makenemo with the -u argument
+   will not invoke wget. For this to be effective the set_local_usp script should be run
+   within the same bash shell as makenemo will be run in. I.e.:
+
+   . ./set_local_usp
+
+   If this is not possible or the makenemo command is run at a later date, then the
+   def_wget script should be run in the bash shell before any attempt to use the -u
+   option of makenemo. I.e.:
+
+   . ./def_wget
+
+   [ Tip: check that the definition/re-definition of wget has worked correctly by
+     issuing this command:
+
+     wget a b c
+
+     If you get messages such as:
+
+       wget a b c
+       --2016-10-27 10:14:52--  http://a/
+       Resolving a... failed: Name or service not known.
+       wget: unable to resolve host address `a'
+       .
+       .
+
+     then the re-definition of wget has not been retained by the current shell. Run
+     . ./def_wget and try again. You should get a response such as:
+
+     wget a b c
+     Expected wget usage: wget src -O dest
+     -O not found. No action taken
+
+     if the environment is set correctly
+   ]
+
+5. The setup is now complete and moving to the CONFIG directory and issuing a command
+   such as:
+
+    ./makenemo -n MYISOMIP -u ISOMIP -m target_arch
+
+   should correctly configure from local archives only.   
+
+   [ Tip: the -u argument ( and hence the need to switch wget to local copies on isolated
+     systems ) is only required when compiling for the first time. Subsequent compilations
+     following code changes in the MYISOMIP directories can be carried out in your
+     default nemo environment. I.e.:
+
+     ./makenemo -n MYISOMIP -m target_arch
+
+     will not require any of the preceeding steps in future sessions.
+   ]
+
Index: branches/2016/dev_CNRS_2016/NEMOGCM/TOOLS/MISCELLANEOUS/make_usp_tar.sh
===================================================================
--- branches/2016/dev_CNRS_2016/NEMOGCM/TOOLS/MISCELLANEOUS/make_usp_tar.sh	(revision 7282)
+++ branches/2016/dev_CNRS_2016/NEMOGCM/TOOLS/MISCELLANEOUS/make_usp_tar.sh	(revision 7282)
@@ -0,0 +1,235 @@
+#!/bin/bash
+#set -x
+set -o posix
+#set -u
+#set -e
+#+
+# ===============
+# make_usp_tar.sh
+# ===============
+# ---------------
+# Fetch and tar a unsupported configuration setup.
+# This script is only needed when target systems do
+# not have wget access to the internet. To configure
+# unsupported conigurations on these systems it will
+# first be necessary to run this script on a system that
+# does have access. Then copy access the resulting tar file,
+# unpack and run the enclosed set_local_uspcfg script in a 
+# bash shell to complete the process. That script redefines/
+# defines wget as a bash function to perform local copies 
+# from this unpacked tarball. If you wish to create a new 
+# configuration based on this local copy of an unsupported 
+# configuration in future sessions then you may need to 
+# redefine wget again before running makenemo with the 
+# appropriate -u setting. A simple script: def_wget is included
+# for such situations.
+# ---------------
+# SYNOPSIS
+# ========
+#  $ make_usp_tar.sh uspcfg.txt target_conf target_dir
+# Note target_dir.tar will be created
+#
+# DESCRIPTION
+# ===========
+# - Extract target configuration details from uspcfg.txt
+# - Create target directory
+# - Recursively use wget to retrieve remote configuration files
+#   into target directory
+# - Copy uspcfg.txt into target directory and alter remote paths
+#   to local (relative) versions
+# - tar target_directory and remove originals
+#----------------------------------------------------------------
+#
+#----------------------------------------------------------------
+# Check the correct number of arguments have been provided
+#----------------------------------------------------------------
+#
+   if [ "$#" != "3" ]; then
+    echo "Expected usage: make_usp_tar.sh uspcfg.txt target_conf target_dir"
+    exit
+   fi
+#
+#----------------------------------------------------------------
+# Check the named uspcfg.txt file exists
+#----------------------------------------------------------------
+#
+   if [ ! -f $1 ]; then
+    echo "named uspcfg.txt file does not exist ($1); attempt abandoned"
+    exit
+   fi
+#
+#----------------------------------------------------------------
+# Check the requested configuration is listed in the named uspcfg.txt file
+#----------------------------------------------------------------
+#
+   inthere=$( grep -c "$2" $1 )
+   if [ "$inthere" -lt 1 ]; then
+    echo "requested configuration is not in named uspcfg.txt file ($2); attempt abandoned"
+    exit
+   fi
+#
+#----------------------------------------------------------------
+# Create the target directory if it does not already exist
+# and cd into it
+#----------------------------------------------------------------
+#
+   if [ ! -d $3 ]; then
+    mkdir $3
+   else
+    echo "target directory already exists; attempt abandoned"
+    exit
+   fi
+   basedir=$(pwd)
+#
+   cd ${3}
+#
+#----------------------------------------------------------------
+# Copy named uspcfg.txt file into target directory
+#----------------------------------------------------------------
+#
+   cp $basedir/$1 .
+#
+#----------------------------------------------------------------
+# Extract information on target configuration and
+# retrieve full file list from remote server
+#----------------------------------------------------------------
+#
+   grep "$2 " $1 > ./cfg.tmp
+#
+   LOCAL_REF=$(cat cfg.tmp | awk 'BEGIN {FS = "#" }{print $2}')
+   TAB=$(cat cfg.tmp | awk 'BEGIN {FS = "#" }{print $3}')
+   REMOTE_CTL=$(cat cfg.tmp | awk 'BEGIN {FS = "#" }{print $4}')
+   wget ${REMOTE_CTL} -O remote_file.list
+#
+#----------------------------------------------------------------
+# Retrieve each remote file and create local directory structure
+# At the same time prepare a local version of the control file
+# by replacing http links with a string that will be replaced later
+# with a local directory path
+#----------------------------------------------------------------
+#
+   if [ -f remote_file.list ] ; then
+    cat remote_file.list | grep -v '^#' |
+     while
+      read remfile locfile
+      do
+       if [ $remfile == 'create_directory' ] ;then
+         mkdir $locfile
+         echo $remfile "                 " $locfile >> local_template.ctl
+       else
+         wget $remfile -O $locfile
+         remfile=$(echo $remfile | sed -e "s;.*$locfile;SET_LOCAL_DIRECTORY_HERE/$locfile;")
+         echo $remfile "  " $locfile >> local_template.ctl
+       fi
+      done
+   else
+    echo "Unable to find remote_file.list. Attempt abandoned"
+    echo "Files may be left in "$(pwd)
+    exit
+   fi
+#
+#----------------------------------------------------------------
+# Construct a modified version of the named uspcfg.txt file
+# First copy across all the other untouched configuration listed
+#----------------------------------------------------------------
+#
+   grep -v "$2 " $1 >  uspcfg_local_template.txt
+#
+#----------------------------------------------------------------
+# Now append the modified entry, replacing http links as before
+#----------------------------------------------------------------
+#
+   cat cfg.tmp | sed -e "s;http.*$;SET_LOCAL_DIRECTORY_HERE/local.ctl;" >> uspcfg_local_template.txt
+#
+#
+#----------------------------------------------------------------
+# Construct a script that can be used later to complete a local installation
+#----------------------------------------------------------------
+#
+cat > set_local_uspcfg << EOF
+#!/bin/bash
+   if [ "\$( echo $SHELL | grep -c bash )" -lt 1 ]; then
+    echo "WARNING: This is only going to be effective in a bash shell"
+    echo "since it redefines the wget command as a bash function."
+    echo "(ignore this comment if you are in a bash shell)"
+   fi
+   basedir=\$(pwd)
+#
+   echo "Enter full path to the CONFIG directory on your target system: "
+   read confdir
+#
+# Edit the local.ctl file to set the local directory path
+#
+   sed -e "s;SET_LOCAL_DIRECTORY_HERE;\$basedir;" local_template.ctl > local.ctl
+#
+# Edit the uspcfg_local.txt file to set the local directory path
+#
+   sed -e "s;SET_LOCAL_DIRECTORY_HERE;\$basedir;" uspcfg_local_template.txt > uspcfg_local.txt
+#
+# Install local versions in the named CONFIG directory
+#
+   if [ -f \$confdir/uspcfg.txt ] && [ ! -L \$confdir/uspcfg.txt ]; then
+    mv \$confdir/uspcfg.txt \$confdir/uspcfg_remote.txt
+    echo "\$confdir/uspcfg.txt moved to \$confdir/uspcfg_remote.txt"
+   fi
+   if [ -f \$confdir/uspcfg_local.txt ]; then
+    echo "Existing uspcfg_local.txt file found in \$confdir"
+    echo "This has been moved to: "\$confdir/uspcfg_local.txt\$\$
+    mv \$confdir/uspcfg_local.txt \$confdir/uspcfg_local.txt\$\$
+   fi
+   mv uspcfg_local.txt \$confdir/uspcfg_local.txt
+   ln -s \$confdir/uspcfg_local.txt \$confdir/uspcfg.txt
+#
+# define/redefine the wget command
+#
+function wget {
+   if [ "\$2" != "-O" ]; then
+    echo "Expected wget usage: wget src -O dest"
+    echo "-O not found. No action taken"
+   else
+    cp \$1 \$3
+   fi
+}
+export -f wget
+EOF
+#----------------------------------------------------------------
+# Construct a script that can be used later to just redefine wget 
+# in bash shell sessions
+#----------------------------------------------------------------
+#
+cat > def_wget << EOFC
+#!/bin/bash
+   if [ "\$( echo $SHELL | grep -c bash )" -lt 1 ]; then
+    echo "WARNING: This is only going to be effective in a bash shell"
+    echo "since it redefines the wget command as a bash function."
+    echo "(ignore this comment if you are in a bash shell)"
+   fi
+function wget {
+   if [ "\$2" != "-O" ]; then
+    echo "Expected wget usage: wget src -O dest"
+    echo "-O not found. No action taken"
+   else
+    cp \$1 \$3
+   fi
+}
+export -f wget
+EOFC
+#----------------------------------------------------------------
+# Make sure these scripts have execute permission
+#----------------------------------------------------------------
+   chmod 755 set_local_uspcfg
+   chmod 755 def_wget
+#
+#----------------------------------------------------------------
+# Tidy up and tar the contents of the downloaded configuration
+#----------------------------------------------------------------
+   rm cfg.tmp
+   cd $basedir
+   tar cvf ${3}.tar $3
+   echo ${3}.tar " file successfully created and prepared for local references. Move this"
+   echo "tar file to your target system, unpack and run the set_local_uspcfg script in a "
+   echo "bash shell to complete the process. This script redefines/defines wget as a bash"
+   echo "function to perform local copies from this unpacked tarball. If you wish to create"
+   echo "a new configuration based on this local copy of an unsupported configuration in "
+   echo "future sessions then you may need to rerun the def_wget script"
+exit
Index: branches/2016/dev_CNRS_2016/NEMOGCM/TOOLS/NESTING/agulhas
===================================================================
--- branches/2016/dev_CNRS_2016/NEMOGCM/TOOLS/NESTING/agulhas	(revision 7280)
+++ branches/2016/dev_CNRS_2016/NEMOGCM/TOOLS/NESTING/agulhas	(revision 7282)
@@ -41,5 +41,5 @@
     N  = 31
     ldbletanh   = .FALSE.
-    ppa2        = 0.0
+    pa2         = 0.0
     ppkth2      = 0.0
     ppacr2      = 0.0
Index: branches/2016/dev_CNRS_2016/NEMOGCM/TOOLS/NESTING/src/agrif_types.f90
===================================================================
--- branches/2016/dev_CNRS_2016/NEMOGCM/TOOLS/NESTING/src/agrif_types.f90	(revision 7280)
+++ branches/2016/dev_CNRS_2016/NEMOGCM/TOOLS/NESTING/src/agrif_types.f90	(revision 7282)
@@ -76,5 +76,5 @@
   NAMELIST /nesting/imin,imax,jmin,jmax,rho,rhot,bathy_update,updated_parent_file      
   !
-  NAMELIST /vertical_grid/ppkth,ppacr,ppdzmin,pphmax,psur,pa0,pa1,N,ldbletanh,ppa2,ppkth2,ppacr2
+  NAMELIST /vertical_grid/ppkth,ppacr,ppdzmin,pphmax,psur,pa0,pa1,N,ldbletanh,pa2,ppkth2,ppacr2
   ! 
   NAMELIST /partial_cells/partial_steps,parent_bathy_meter,parent_batmet_name,e3zps_min,e3zps_rat      
