Index: /NEMO/branches/2019/dev_r11943_MERGE_2019/cfgs/SHARED/field_def_nemo-oce.xml
===================================================================
--- /NEMO/branches/2019/dev_r11943_MERGE_2019/cfgs/SHARED/field_def_nemo-oce.xml	(revision 12371)
+++ /NEMO/branches/2019/dev_r11943_MERGE_2019/cfgs/SHARED/field_def_nemo-oce.xml	(revision 12372)
@@ -388,4 +388,5 @@
           <field id="hflx_rain_cea" long_name="heat flux due to rainfall"                                standard_name="temperature_flux_due_to_rainfall_expressed_as_heat_flux_into_sea_water"        unit="W/m2"     />
           <field id="hflx_evap_cea" long_name="heat flux due to evaporation"                             standard_name="temperature_flux_due_to_evaporation_expressed_as_heat_flux_out_of_sea_water"   unit="W/m2"     />
+          <field id="hflx_prec_cea" long_name="heat flux due to all precip"                              standard_name="temperature_flux_due_to_all_precip_expressed_as_heat_flux_into_sea_water"      unit="W/m2"     />
           <field id="hflx_snow_cea" long_name="heat flux due to snow falling"                            standard_name="heat_flux_onto_ocean_and_ice_due_to_snow_thermodynamics"                       unit="W/m2"     />
           <field id="hflx_snow_ai_cea" long_name="heat flux due to snow falling over ice"                standard_name="heat_flux_onto_ice_due_to_snow_thermodynamics"                                 unit="W/m2"     />
Index: /NEMO/branches/2019/dev_r11943_MERGE_2019/cfgs/SHARED/grid_def_nemo.xml
===================================================================
--- /NEMO/branches/2019/dev_r11943_MERGE_2019/cfgs/SHARED/grid_def_nemo.xml	(revision 12371)
+++ /NEMO/branches/2019/dev_r11943_MERGE_2019/cfgs/SHARED/grid_def_nemo.xml	(revision 12372)
@@ -102,4 +102,28 @@
          <axis axis_ref="depthw"  />
          <axis axis_ref="basin" />
+       </grid>
+
+       <grid id="grid_ptr_W_GLO">
+         <domain  domain_ref="ptr" />
+         <axis axis_ref="depthw"  />
+         <scalar>
+           <extract_axis position="0" />
+         </scalar>
+       </grid>
+
+       <grid id="grid_ptr_W_ATL">
+         <domain  domain_ref="ptr" />
+         <axis axis_ref="depthw"  />
+         <scalar>
+           <extract_axis position="1" />
+         </scalar>
+       </grid>
+
+       <grid id="grid_ptr_W_IND">
+         <domain  domain_ref="ptr" />
+         <axis axis_ref="depthw"  />
+         <scalar>
+           <extract_axis position="2" />
+         </scalar>
        </grid>
 
Index: /NEMO/branches/2019/dev_r11943_MERGE_2019/src/OCE/DYN/dynatf.F90
===================================================================
--- /NEMO/branches/2019/dev_r11943_MERGE_2019/src/OCE/DYN/dynatf.F90	(revision 12371)
+++ /NEMO/branches/2019/dev_r11943_MERGE_2019/src/OCE/DYN/dynatf.F90	(revision 12372)
@@ -102,5 +102,5 @@
       REAL(wp) ::   zue3a, zue3n, zue3b, zcoef    ! local scalars
       REAL(wp) ::   zve3a, zve3n, zve3b, z1_2dt   !   -      -
-      REAL(wp), ALLOCATABLE, DIMENSION(:,:)   ::   zue, zve
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:)   ::   zue, zve, zwfld
       REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   ze3t_f, ze3u_f, ze3v_f, zua, zva 
       !!----------------------------------------------------------------------
@@ -191,26 +191,17 @@
             ! Time-filtered scale factor at t-points
             ! ----------------------------------------------------
-            ALLOCATE( ze3t_f(jpi,jpj,jpk) )
+            ALLOCATE( ze3t_f(jpi,jpj,jpk), zwfld(jpi,jpj) )
             DO jk = 1, jpkm1
                ze3t_f(:,:,jk) = pe3t(:,:,jk,Kmm) + atfp * ( pe3t(:,:,jk,Kbb) - 2._wp * pe3t(:,:,jk,Kmm) + pe3t(:,:,jk,Kaa) )
             END DO
             ! Add volume filter correction: compatibility with tracer advection scheme
-            ! => time filter + conservation correction (only at the first level)
+            ! => time filter + conservation correction
             zcoef = atfp * rdt * r1_rau0
-
-            ze3t_f(:,:,1) = ze3t_f(:,:,1) - zcoef * ( emp_b(:,:) - emp(:,:) ) * tmask(:,:,1)
-
-            IF ( ln_rnf ) THEN
-               IF( ln_rnf_depth ) THEN
-                  DO_3D_11_11( 1, jpkm1 )
-                     IF( jk <=  nk_rnf(ji,jj)  ) THEN
-                         ze3t_f(ji,jj,jk) =   ze3t_f(ji,jj,jk) - zcoef *  ( - rnf_b(ji,jj) + rnf(ji,jj) ) &
-                                &          * ( pe3t(ji,jj,jk,Kmm) / h_rnf(ji,jj) ) * tmask(ji,jj,jk)
-                     ENDIF
-                  END_3D
-               ELSE
-                  ze3t_f(:,:,1) = ze3t_f(:,:,1) - zcoef *  ( -rnf_b(:,:) + rnf(:,:))*tmask(:,:,1)
-               ENDIF
-            END IF
+            zwfld(:,:) = emp_b(:,:) - emp(:,:)
+            IF ( ln_rnf ) zwfld(:,:) =  zwfld(:,:) - ( rnf_b(:,:) - rnf(:,:) )
+            DO jk = 1, jpkm1
+               ze3t_f(:,:,jk) = ze3t_f(:,:,jk) - zcoef * zwfld(:,:) * tmask(:,:,jk) &
+                              &                        * pe3t(:,:,jk,Kmm) / ( ht(:,:) + 1._wp - ssmask(:,:) ) 
+            END DO
             !
             ! ice shelf melting (deal separately as it can be in depth)
@@ -254,5 +245,5 @@
             ENDIF
             !
-            DEALLOCATE( ze3t_f )
+            DEALLOCATE( ze3t_f, zwfld )
          ENDIF
          !
Index: /NEMO/branches/2019/dev_r11943_MERGE_2019/src/OCE/ISF/isfdynatf.F90
===================================================================
--- /NEMO/branches/2019/dev_r11943_MERGE_2019/src/OCE/ISF/isfdynatf.F90	(revision 12371)
+++ /NEMO/branches/2019/dev_r11943_MERGE_2019/src/OCE/ISF/isfdynatf.F90	(revision 12372)
@@ -14,5 +14,5 @@
 
    USE phycst , ONLY: r1_rau0         ! physical constant
-   USE dom_oce, ONLY: e3t, r1_e1e2t   ! time and space domain
+   USE dom_oce, ONLY: tmask, ssmask, ht, e3t, r1_e1e2t   ! time and space domain
 
    USE in_out_manager
@@ -77,14 +77,10 @@
       !
       ! compute fwf conservation correction
-      zfwfinc(:,:) = pcoef * ( pfwf_b(:,:) - pfwf(:,:) ) / phtbl(:,:) * r1_rau0
+      zfwfinc(:,:) = pcoef * ( pfwf_b(:,:) - pfwf(:,:) ) / ( ht(:,:) + 1._wp - ssmask(:,:) ) * r1_rau0
       !
-      ! add the increment in the tbl
-      DO_3D_11_11( 1, jpkm1 )
-         IF( ktop(ji,jj) <= jk .AND. jk < kbot(ji,jj)  ) THEN
-            pe3t_f(ji,jj,jk) = pe3t_f(ji,jj,jk) - zfwfinc(ji,jj) * e3t(ji,jj,jk,Kmm)
-         ELSEIF ( jk == kbot(ji,jj) ) THEN
-            pe3t_f(ji,jj,jk) = pe3t_f(ji,jj,jk) - zfwfinc(ji,jj) * e3t(ji,jj,jk,Kmm) * pfrac(ji,jj)
-         ENDIF
-      END_3D
+      ! add the increment
+      DO jk = 1, jpkm1
+         pe3t_f(:,:,jk) = pe3t_f(:,:,jk) - tmask(:,:,jk) * zfwfinc(:,:) * e3t(:,:,jk,Kmm)
+      END DO
       !
    END SUBROUTINE isf_dynatf_mlt
Index: /NEMO/branches/2019/dev_r11943_MERGE_2019/src/OCE/SBC/fldread.F90
===================================================================
--- /NEMO/branches/2019/dev_r11943_MERGE_2019/src/OCE/SBC/fldread.F90	(revision 12371)
+++ /NEMO/branches/2019/dev_r11943_MERGE_2019/src/OCE/SBC/fldread.F90	(revision 12372)
@@ -571,6 +571,6 @@
 
       REAL(wp), DIMENSION(:,:,:), INTENT(in   ) ::   pdta_read       ! data read in bdy file
-      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   pdta_read_z     ! depth of the data read in bdy file
-      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   pdta_read_dz    ! thickness of the levels in bdy file
+      REAL(wp), DIMENSION(:,:,:), INTENT(in   ) ::   pdta_read_z     ! depth of the data read in bdy file
+      REAL(wp), DIMENSION(:,:,:), INTENT(in   ) ::   pdta_read_dz    ! thickness of the levels in bdy file
       REAL(wp), DIMENSION(:,:,:), INTENT(  out) ::   pdta            ! output field on model grid (2 dimensional)
       REAL(wp)                  , INTENT(in   ) ::   pfv             ! fillvalue of the data read in bdy file
@@ -580,11 +580,11 @@
       INTEGER, OPTIONAL         , INTENT(in   ) ::   Kmm             ! ocean time level index
       !!
-      INTEGER                                   ::   ipi             ! length of boundary data on local process
-      INTEGER                                   ::   ipkb            ! number of vertical levels in boundary data file
-      INTEGER                                   ::   jb, ji, jj, jk, jkb   ! loop counters
-      REAL(wp)                                  ::   zcoef
-      REAL(wp)                                  ::   zl, zi, zh      ! tmp variable for current depth and interpolation factor
-      REAL(wp)                                  ::   zfv_alt, ztrans, ztrans_new ! fillvalue and alternative -ABS(pfv)
-      REAL(wp), DIMENSION(jpk)                  ::   zdepth, zdhalf  ! level and half-level depth
+      INTEGER                  ::   ipi                 ! length of boundary data on local process
+      INTEGER                  ::   ipkb                ! number of vertical levels in boundary data file
+      INTEGER                  ::   ipkmax              ! number of vertical levels in boundary data file where no mask
+      INTEGER                  ::   jb, ji, jj, jk, jkb ! loop counters
+      REAL(wp)                 ::   zcoef, zi           ! 
+      REAL(wp)                 ::   ztrans, ztrans_new  ! transports
+      REAL(wp), DIMENSION(jpk) ::   zdepth, zdhalf      ! level and half-level depth
       !!---------------------------------------------------------------------
      
@@ -592,23 +592,20 @@
       ipkb = SIZE( pdta_read, 3 )
       
-      zfv_alt = -ABS(pfv)  ! set _FillValue < 0 as we make use of MAXVAL and MAXLOC later
-      !
-      WHERE( pdta_read == pfv )
-         pdta_read_z  = zfv_alt ! safety: put fillvalue into external depth field so consistent with data
-         pdta_read_dz = 0._wp   ! safety: put 0._wp into external thickness factors to ensure transport is correct
-      ENDWHERE
-      
       DO jb = 1, ipi
          ji = idx_bdy(kbdy)%nbi(jb,kgrd)
          jj = idx_bdy(kbdy)%nbj(jb,kgrd)
-         zh  = SUM(pdta_read_dz(jb,1,:) )
-         !
+         !
+         ! --- max jk where input data /= FillValue --- !
+         ipkmax = 1
+         DO jkb = 2, ipkb
+            IF( pdta_read(jb,1,jkb) /= pfv )   ipkmax = MAX( ipkmax, jkb )
+         END DO
+         !
+         ! --- calculate depth at t,u,v points --- !
          SELECT CASE( kgrd )                         
-         CASE(1)
-            ! depth of T points:
+         CASE(1)            ! depth of T points:
             zdepth(:) = gdept(ji,jj,:,Kmm)
-         CASE(2)
-            ! depth of U points: we must not use gdept_n as we don't want to do a communication
-            !   --> copy what is done for gdept_n in domvvl...
+         CASE(2)            ! depth of U points: we must not use gdept_n as we don't want to do a communication
+            !                 --> copy what is done for gdept_n in domvvl...
             zdhalf(1) = 0.0_wp
             zdepth(1) = 0.5_wp * e3uw(ji,jj,1,Kmm)
@@ -620,10 +617,9 @@
                zcoef = ( umask(ji,jj,jk) - wumask(ji,jj,jk) )
                zdhalf(jk) = zdhalf(jk-1) + e3u(ji,jj,jk-1,Kmm)
-               zdepth(jk) =      zcoef  * ( zdhalf(jk  ) + 0.5 * e3uw(ji,jj,jk,Kmm))  &
-                  &         + (1-zcoef) * ( zdepth(jk-1) + e3uw(ji,jj,jk,Kmm))
+               zdepth(jk) =          zcoef  * ( zdhalf(jk  ) + 0.5_wp * e3uw(ji,jj,jk,Kmm))  &
+                  &         + (1._wp-zcoef) * ( zdepth(jk-1) +          e3uw(ji,jj,jk,Kmm))
             END DO
-         CASE(3)
-            ! depth of V points: we must not use gdept_n as we don't want to do a communication
-            !   --> copy what is done for gdept_n in domvvl...
+         CASE(3)            ! depth of V points: we must not use gdept_n as we don't want to do a communication
+            !                 --> copy what is done for gdept_n in domvvl...
             zdhalf(1) = 0.0_wp
             zdepth(1) = 0.5_wp * e3vw(ji,jj,1,Kmm)
@@ -635,20 +631,20 @@
                zcoef = ( vmask(ji,jj,jk) - wvmask(ji,jj,jk) )
                zdhalf(jk) = zdhalf(jk-1) + e3v(ji,jj,jk-1,Kmm)
-               zdepth(jk) =      zcoef  * ( zdhalf(jk  ) + 0.5 * e3vw(ji,jj,jk,Kmm))  &
-                  &         + (1-zcoef) * ( zdepth(jk-1) + e3vw(ji,jj,jk,Kmm))
+               zdepth(jk) =          zcoef  * ( zdhalf(jk  ) + 0.5_wp * e3vw(ji,jj,jk,Kmm))  &
+                  &         + (1._wp-zcoef) * ( zdepth(jk-1) +          e3vw(ji,jj,jk,Kmm))
             END DO
          END SELECT
-         !
-         DO jk = 1, jpk                      
-            IF(     zdepth(jk) < pdta_read_z(jb,1,          1) ) THEN                ! above the first level of external data
-               pdta(jb,1,jk) =  pdta_read(jb,1,1)
-            ELSEIF( zdepth(jk) > pdta_read_z(jb,1,ipkb) ) THEN                       ! below the last level of external data 
-               pdta(jb,1,jk) =  pdta_read(jb,1,MAXLOC(pdta_read_z(jb,1,:),1))
-            ELSE                                                             ! inbetween: vertical interpolation between jkb & jkb+1
-               DO jkb = 1, ipkb-1                                            ! when  gdept_n(jkb) < zdepth(jk) < gdept_n(jkb+1)
-                  IF( ( ( zdepth(jk) - pdta_read_z(jb,1,jkb) ) * ( zdepth(jk) - pdta_read_z(jb,1,jkb+1) ) <= 0._wp ) &
-                     &    .AND. ( pdta_read_z(jb,1,jkb+1) /= zfv_alt) ) THEN   ! linear interpolation between 2 levels
+         !         
+         ! --- interpolate bdy data on the model grid --- !
+         DO jk = 1, jpk
+            IF(     zdepth(jk) <= pdta_read_z(jb,1,1)      ) THEN   ! above the first level of external data
+               pdta(jb,1,jk) = pdta_read(jb,1,1)
+            ELSEIF( zdepth(jk) >  pdta_read_z(jb,1,ipkmax) ) THEN   ! below the last level of external data /= FillValue
+               pdta(jb,1,jk) = pdta_read(jb,1,ipkmax)
+            ELSE                                                    ! inbetween: vertical interpolation between jkb & jkb+1
+               DO jkb = 1, ipkmax-1
+                  IF( ( ( zdepth(jk) - pdta_read_z(jb,1,jkb) ) * ( zdepth(jk) - pdta_read_z(jb,1,jkb+1) ) ) <= 0._wp ) THEN ! linear interpolation between 2 levels
                      zi = ( zdepth(jk) - pdta_read_z(jb,1,jkb) ) / ( pdta_read_z(jb,1,jkb+1) - pdta_read_z(jb,1,jkb) )
-                     pdta(jb,1,jk) = pdta_read(jb,1,jkb) + ( pdta_read  (jb,1,jkb+1) - pdta_read  (jb,1,jkb) ) * zi
+                     pdta(jb,1,jk) = pdta_read(jb,1,jkb) + zi * ( pdta_read(jb,1,jkb+1) - pdta_read(jb,1,jkb) )
                   ENDIF
                END DO
@@ -657,15 +653,26 @@
          !
       END DO   ! ipi
-      
-      IF(kgrd == 2) THEN                                  ! do we need to adjust the transport term?
+
+      ! --- mask data and adjust transport --- !
+      SELECT CASE( kgrd )                         
+
+      CASE(1)                                 ! mask data (probably unecessary)
          DO jb = 1, ipi
             ji = idx_bdy(kbdy)%nbi(jb,kgrd)
             jj = idx_bdy(kbdy)%nbj(jb,kgrd)
-            zh  = SUM(pdta_read_dz(jb,1,:) )
+            DO jk = 1, jpk                      
+               pdta(jb,1,jk) = pdta(jb,1,jk) * tmask(ji,jj,jk)
+            END DO
+         END DO
+         
+      CASE(2)                                 ! adjust the U-transport term
+         DO jb = 1, ipi
+            ji = idx_bdy(kbdy)%nbi(jb,kgrd)
+            jj = idx_bdy(kbdy)%nbj(jb,kgrd)
             ztrans = 0._wp
+            DO jkb = 1, ipkb                              ! calculate transport on input grid
+               IF( pdta_read(jb,1,jkb) /= pfv )   ztrans = ztrans + pdta_read(jb,1,jkb) * pdta_read_dz(jb,1,jkb)
+            ENDDO
             ztrans_new = 0._wp
-            DO jkb = 1, ipkb                              ! calculate transport on input grid
-               ztrans     = ztrans     + pdta_read(jb,1,jkb) * pdta_read_dz(jb, 1,jkb)
-            ENDDO
             DO jk = 1, jpk                                ! calculate transport on model grid
                ztrans_new = ztrans_new +      pdta(jb,1,jk ) * e3u(ji,jj,jk,Kmm ) * umask(ji,jj,jk)
@@ -674,33 +681,31 @@
                IF(ldtotvel) THEN ! bdy data are total velocity so adjust bt transport term to match input data
                   pdta(jb,1,jk) = ( pdta(jb,1,jk) + ( ztrans - ztrans_new ) * r1_hu(ji,jj,Kmm) ) * umask(ji,jj,jk)
-               ELSE ! we're just dealing with bc velocity so bt transport term should sum to zero
+               ELSE              ! we're just dealing with bc velocity so bt transport term should sum to zero
                   pdta(jb,1,jk) =   pdta(jb,1,jk) + (  0._wp - ztrans_new ) * r1_hu(ji,jj,Kmm)   * umask(ji,jj,jk)
                ENDIF
             ENDDO
          ENDDO
-      ENDIF
-      
-      IF(kgrd == 3) THEN                                  ! do we need to adjust the transport term?
+
+      CASE(3)                                 ! adjust the V-transport term
          DO jb = 1, ipi
             ji = idx_bdy(kbdy)%nbi(jb,kgrd)
             jj = idx_bdy(kbdy)%nbj(jb,kgrd)
-            zh  = SUM(pdta_read_dz(jb,1,:) )
             ztrans = 0._wp
+            DO jkb = 1, ipkb                              ! calculate transport on input grid
+               IF( pdta_read(jb,1,jkb) /= pfv )   ztrans = ztrans + pdta_read(jb,1,jkb) * pdta_read_dz(jb,1,jkb)
+            ENDDO
             ztrans_new = 0._wp
-            DO jkb = 1, ipkb                              ! calculate transport on input grid
-               ztrans     = ztrans     + pdta_read(jb,1,jkb) * pdta_read_dz(jb, 1,jkb)
-            ENDDO
             DO jk = 1, jpk                                ! calculate transport on model grid
-               ztrans_new = ztrans_new +      pdta(jb,1,jk ) *        e3v(ji,jj,jk,Kmm ) * vmask(ji,jj,jk)
+               ztrans_new = ztrans_new +      pdta(jb,1,jk ) * e3v(ji,jj,jk,Kmm ) * vmask(ji,jj,jk)
             ENDDO
             DO jk = 1, jpk                                ! make transport correction
                IF(ldtotvel) THEN ! bdy data are total velocity so adjust bt transport term to match input data
                   pdta(jb,1,jk) = ( pdta(jb,1,jk) + ( ztrans - ztrans_new ) * r1_hv(ji,jj,Kmm) ) * vmask(ji,jj,jk)
-               ELSE ! we're just dealing with bc velocity so bt transport term should sum to zero
+               ELSE              ! we're just dealing with bc velocity so bt transport term should sum to zero
                   pdta(jb,1,jk) =   pdta(jb,1,jk) + (  0._wp - ztrans_new ) * r1_hv(ji,jj,Kmm)   * vmask(ji,jj,jk)
                ENDIF
             ENDDO
          ENDDO
-      ENDIF
+      END SELECT
       
    END SUBROUTINE fld_bdy_interp
Index: /NEMO/branches/2019/dev_r11943_MERGE_2019/src/OCE/TRA/traatf.F90
===================================================================
--- /NEMO/branches/2019/dev_r11943_MERGE_2019/src/OCE/TRA/traatf.F90	(revision 12371)
+++ /NEMO/branches/2019/dev_r11943_MERGE_2019/src/OCE/TRA/traatf.F90	(revision 12372)
@@ -252,5 +252,5 @@
       INTEGER  ::   ji, jj, jk, jn              ! dummy loop indices
       REAL(wp) ::   zfact, zfact1, ztc_a , ztc_n , ztc_b , ztc_f , ztc_d    ! local scalar
-      REAL(wp) ::   zfact2, ze3t_b, ze3t_n, ze3t_a, ze3t_f, ze3t_d   !   -      -
+      REAL(wp) ::   zfact2, ze3t_b, ze3t_n, ze3t_a, ze3t_f, ze3t_d, zscale  !   -      -
       REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:) ::   ztrd_atf
       !!----------------------------------------------------------------------
@@ -295,20 +295,16 @@
             ztc_f  = ztc_n  + atfp * ztc_d
             !
+            ! Add asselin correction on scale factors:
+            zscale = tmask(ji,jj,jk) * e3t(ji,jj,jk,Kmm) / ( ht(ji,jj) + 1._wp - ssmask(ji,jj) ) 
+            ze3t_f = ze3t_f - zfact2 * zscale * ( emp_b(ji,jj) - emp(ji,jj) ) 
+            IF ( ll_rnf ) ze3t_f = ze3t_f + zfact2 * zscale * (    rnf_b(ji,jj) -    rnf(ji,jj) ) 
+            IF ( ll_isf ) THEN
+               IF ( ln_isfcav_mlt ) ze3t_f = ze3t_f - zfact2 * zscale * ( fwfisf_cav_b(ji,jj) - fwfisf_cav(ji,jj) )
+               IF ( ln_isfpar_mlt ) ze3t_f = ze3t_f - zfact2 * zscale * ( fwfisf_par_b(ji,jj) - fwfisf_par(ji,jj) )
+            ENDIF
+            !
             IF( jk == mikt(ji,jj) ) THEN           ! first level 
-               ze3t_f = ze3t_f - zfact2 * ( emp_b(ji,jj)    - emp(ji,jj) )
                ztc_f  = ztc_f  - zfact1 * ( psbc_tc(ji,jj,jn) - psbc_tc_b(ji,jj,jn) )
             ENDIF
-            IF( ln_rnf_depth ) THEN
-               ! Rivers are not just at the surface must go down to nk_rnf(ji,jj)
-               IF( jk <= nk_rnf(ji,jj)  ) THEN
-                  ze3t_f = ze3t_f - zfact2 * ( - (rnf_b(ji,jj) - rnf(ji,jj)   )  ) &
-              &                            * ( e3t(ji,jj,jk,Kmm) / h_rnf(ji,jj) ) 
-               ENDIF
-            ELSE
-               IF( jk == 1 ) THEN           ! first level 
-                  ze3t_f = ze3t_f - zfact2 * ( - (rnf_b(ji,jj)    - rnf(ji,jj)   ) ) 
-               ENDIF
-            ENDIF
-            !
             !
             ! solar penetration (temperature only)
@@ -331,6 +327,4 @@
                      ztc_f  = ztc_f  - zfact1 * ( risf_cav_tsc(ji,jj,jn) - risf_cav_tsc_b(ji,jj,jn) ) &
                         &                     * e3t(ji,jj,jk,Kmm) / rhisf_tbl_cav(ji,jj)
-                     ze3t_f = ze3t_f - zfact2 * ( fwfisf_cav_b(ji,jj) - fwfisf_cav(ji,jj) )           &
-                        &                     * e3t(ji,jj,jk,Kmm) / rhisf_tbl_cav(ji,jj)
                   END IF
                   ! level partially include in Losch_2008 ice shelf boundary layer 
@@ -338,6 +332,4 @@
                      ztc_f  = ztc_f  - zfact1 * ( risf_cav_tsc(ji,jj,jn) - risf_cav_tsc_b(ji,jj,jn) )  &
                             &                 * e3t(ji,jj,jk,Kmm) / rhisf_tbl_cav(ji,jj) * rfrac_tbl_cav(ji,jj)
-                     ze3t_f = ze3t_f - zfact2 * ( fwfisf_cav_b(ji,jj) - fwfisf_cav(ji,jj) )            &
-                        &                     * e3t(ji,jj,jk,Kmm) / rhisf_tbl_cav(ji,jj) * rfrac_tbl_cav(ji,jj)
                   END IF
                END IF
@@ -349,6 +341,4 @@
                      ztc_f  = ztc_f  - zfact1 * ( risf_par_tsc(ji,jj,jn) - risf_par_tsc_b(ji,jj,jn) )  &
                             &                 * e3t(ji,jj,jk,Kmm) / rhisf_tbl_par(ji,jj)
-                     ze3t_f = ze3t_f - zfact2 * ( fwfisf_par_b(ji,jj) - fwfisf_par(ji,jj) )            &
-                        &                     * e3t(ji,jj,jk,Kmm) / rhisf_tbl_par(ji,jj)
                   END IF
                   ! level partially include in Losch_2008 ice shelf boundary layer 
@@ -356,6 +346,4 @@
                      ztc_f  = ztc_f  - zfact1 * ( risf_par_tsc(ji,jj,jn) - risf_par_tsc_b(ji,jj,jn) )  &
                             &                 * e3t(ji,jj,jk,Kmm) / rhisf_tbl_par(ji,jj) * rfrac_tbl_par(ji,jj)
-                     ze3t_f = ze3t_f - zfact2 * ( fwfisf_par_b(ji,jj) - fwfisf_par(ji,jj) )            &
-                        &                     * e3t(ji,jj,jk,Kmm) / rhisf_tbl_par(ji,jj) * rfrac_tbl_par(ji,jj)
                   END IF
                END IF
@@ -367,4 +355,5 @@
                   IF ( ln_rstart .AND. kt == nit000+1 ) THEN
                      ztc_f  = ztc_f  + zfact1 * risfcpl_tsc(ji,jj,jk,jn) * r1_e1e2t(ji,jj)
+                     ! Shouldn't volume increment be spread according thanks to zscale  ?
                      ze3t_f = ze3t_f - zfact1 * risfcpl_vol(ji,jj,jk   ) * r1_e1e2t(ji,jj)
                   END IF
Index: /NEMO/branches/2019/dev_r11943_MERGE_2019/src/TOP/PISCES/P4Z/p4zprod.F90
===================================================================
--- /NEMO/branches/2019/dev_r11943_MERGE_2019/src/TOP/PISCES/P4Z/p4zprod.F90	(revision 12371)
+++ /NEMO/branches/2019/dev_r11943_MERGE_2019/src/TOP/PISCES/P4Z/p4zprod.F90	(revision 12372)
@@ -323,5 +323,5 @@
        CALL iom_put( "LDlight" , zprdia (:,:,:) / (zprmaxd(:,:,:) + rtrn) * tmask(:,:,:)   )
        CALL iom_put( "TPP"     , ( zprorcan(:,:,:) + zprorcad(:,:,:) ) * zfact * tmask(:,:,:)  )  ! total primary production
-       CALL iom_put( "TPNEW"   , ( zprorcan(:,:,:) + zprorcad(:,:,:) ) * zfact * tmask(:,:,:)  ) ! total new production
+       CALL iom_put( "TPNEW"   , ( zpronewn(:,:,:) + zpronewd(:,:,:) ) * zfact * tmask(:,:,:)  ) ! total new production
        CALL iom_put( "TPBFE"   , ( zprofen(:,:,:) + zprofed(:,:,:) ) * zfact * tmask(:,:,:)  )  ! total biogenic iron production
        CALL iom_put( "tintpp"  , tpp * zfact )  !  global total integrated primary production molC/s
