Index: branches/2017/dev_merge_2017/NEMOGCM/NEMO/NST_SRC/agrif_lim3_interp.F90
===================================================================
--- branches/2017/dev_merge_2017/NEMOGCM/NEMO/NST_SRC/agrif_lim3_interp.F90	(revision 9167)
+++ branches/2017/dev_merge_2017/NEMOGCM/NEMO/NST_SRC/agrif_lim3_interp.F90	(revision 9454)
@@ -24,4 +24,5 @@
    USE ice
    USE agrif_ice
+   USE phycst , ONLY: rt0
    
    IMPLICIT NONE
@@ -183,5 +184,5 @@
       ELSE               ! child grid
          !
-         IF( nbghostcells > 1 ) THEN   ! ==> The easiest interpolation is used
+!         IF( nbghostcells > 1 ) THEN   ! ==> The easiest interpolation is used
             !
             jm = 1
@@ -214,144 +215,146 @@
             END DO
             !
-         ELSE                          ! ==> complex interpolation (only one ghost cell available)
-            !! Use a more complex interpolation since we mix solutions over a couple of grid points
-            !! it is advised to use it for fields modified by high order schemes (e.g. advection UM5...)
-            !! clem: for some reason (I don't know why), the following lines do not work 
-            !!       with mpp (or in realistic configurations?). It makes the model crash
-            !        I think there is an issue with Agrif_SpecialValue here (not taken into account properly)
-            ! record ztab
-            jm = 1
-            DO jl = 1, jpl
-               ztab(:,:,jm  ) = a_i (:,:,jl)
-               ztab(:,:,jm+1) = v_i (:,:,jl)
-               ztab(:,:,jm+2) = v_s (:,:,jl)
-               ztab(:,:,jm+3) = sv_i(:,:,jl)
-               ztab(:,:,jm+4) = oa_i(:,:,jl)
-               ztab(:,:,jm+5) = a_ip(:,:,jl)
-               ztab(:,:,jm+6) = v_ip(:,:,jl)
-               ztab(:,:,jm+7) = t_su(:,:,jl)
-               jm = jm + 8
-               DO jk = 1, nlay_s
-                  ztab(:,:,jm) = e_s(:,:,jk,jl)
-                  jm = jm + 1
-               END DO
-               DO jk = 1, nlay_i
-                  ztab(:,:,jm) = e_i(:,:,jk,jl)
-                  jm = jm + 1
-               END DO
-               !
-            END DO
-            !
-            ! borders of the domain
-            western_side  = (nb == 1).AND.(ndir == 1)  ;  eastern_side  = (nb == 1).AND.(ndir == 2)
-            southern_side = (nb == 2).AND.(ndir == 1)  ;  northern_side = (nb == 2).AND.(ndir == 2)
-            !
-            ! spatial smoothing
-            zrhox = Agrif_Rhox()
-            z1 =      ( zrhox - 1. ) * 0.5 
-            z3 =      ( zrhox - 1. ) / ( zrhox + 1. )
-            z6 = 2. * ( zrhox - 1. ) / ( zrhox + 1. )
-            z7 =    - ( zrhox - 1. ) / ( zrhox + 3. )
-            z2 = 1. - z1
-            z4 = 1. - z3
-            z5 = 1. - z6 - z7
-            !
-            ! Remove corners
-            imin = i1  ;  imax = i2  ;  jmin = j1  ;  jmax = j2
-            IF( (nbondj == -1) .OR. (nbondj == 2) )   jmin = 3
-            IF( (nbondj == +1) .OR. (nbondj == 2) )   jmax = nlcj-2
-            IF( (nbondi == -1) .OR. (nbondi == 2) )   imin = 3
-            IF( (nbondi == +1) .OR. (nbondi == 2) )   imax = nlci-2
-
-            ! smoothed fields
-            IF( eastern_side ) THEN
-               ztab(nlci,j1:j2,:) = z1 * ptab(nlci,j1:j2,:) + z2 * ptab(nlci-1,j1:j2,:)
-               DO jj = jmin, jmax
-                  rswitch = 0.
-                  IF( u_ice(nlci-2,jj) > 0._wp ) rswitch = 1.
-                  ztab(nlci-1,jj,:) = ( 1. - umask(nlci-2,jj,1) ) * ztab(nlci,jj,:)  &
-                     &                +      umask(nlci-2,jj,1)   *  &
-                     &                ( ( 1. - rswitch ) * ( z4 * ztab(nlci,jj,:)   + z3 * ztab(nlci-2,jj,:) )  &
-                     &                  +      rswitch   * ( z6 * ztab(nlci-2,jj,:) + z5 * ztab(nlci,jj,:) + z7 * ztab(nlci-3,jj,:) ) )
-                  ztab(nlci-1,jj,:) = ztab(nlci-1,jj,:) * tmask(nlci-1,jj,1)
-               END DO
-            ENDIF
-            ! 
-            IF( northern_side ) THEN
-               ztab(i1:i2,nlcj,:) = z1 * ptab(i1:i2,nlcj,:) + z2 * ptab(i1:i2,nlcj-1,:)
-               DO ji = imin, imax
-                  rswitch = 0.
-                  IF( v_ice(ji,nlcj-2) > 0._wp ) rswitch = 1.
-                  ztab(ji,nlcj-1,:) = ( 1. - vmask(ji,nlcj-2,1) ) * ztab(ji,nlcj,:)  &
-                     &                +      vmask(ji,nlcj-2,1)   *  &
-                     &                ( ( 1. - rswitch ) * ( z4 * ztab(ji,nlcj,:)   + z3 * ztab(ji,nlcj-2,:) ) &
-                     &                  +      rswitch   * ( z6 * ztab(ji,nlcj-2,:) + z5 * ztab(ji,nlcj,:) + z7 * ztab(ji,nlcj-3,:) ) )
-                  ztab(ji,nlcj-1,:) = ztab(ji,nlcj-1,:) * tmask(ji,nlcj-1,1)
-               END DO
-            END IF
-            !
-            IF( western_side) THEN
-               ztab(1,j1:j2,:) = z1 * ptab(1,j1:j2,:) + z2 * ptab(2,j1:j2,:)
-               DO jj = jmin, jmax
-                  rswitch = 0.
-                  IF( u_ice(2,jj) < 0._wp ) rswitch = 1.
-                  ztab(2,jj,:) = ( 1. - umask(2,jj,1) ) * ztab(1,jj,:)  &
-                     &           +      umask(2,jj,1)   *   &
-                     &           ( ( 1. - rswitch ) * ( z4 * ztab(1,jj,:) + z3 * ztab(3,jj,:) ) &
-                     &             +      rswitch   * ( z6 * ztab(3,jj,:) + z5 * ztab(1,jj,:) + z7 * ztab(4,jj,:) ) )
-                  ztab(2,jj,:) = ztab(2,jj,:) * tmask(2,jj,1)
-               END DO
-            ENDIF
-            !
-            IF( southern_side ) THEN
-               ztab(i1:i2,1,:) = z1 * ptab(i1:i2,1,:) + z2 * ptab(i1:i2,2,:)
-               DO ji = imin, imax
-                  rswitch = 0.
-                  IF( v_ice(ji,2) < 0._wp ) rswitch = 1.
-                  ztab(ji,2,:) = ( 1. - vmask(ji,2,1) ) * ztab(ji,1,:)  &
-                     &           +      vmask(ji,2,1)   *  &
-                     &           ( ( 1. - rswitch ) * ( z4 * ztab(ji,1,:) + z3 * ztab(ji,3,:) ) &
-                     &             +      rswitch   * ( z6 * ztab(ji,3,:) + z5 * ztab(ji,1,:) + z7 * ztab(ji,4,:) ) )
-                  ztab(ji,2,:) = ztab(ji,2,:) * tmask(ji,2,1)
-               END DO
-            END IF
-            !
-            ! Treatment of corners
-            IF( (eastern_side) .AND. ((nbondj == -1).OR.(nbondj == 2)) )  ztab(nlci-1,2,:)      = ptab(nlci-1,2,:)      ! East south
-            IF( (eastern_side) .AND. ((nbondj ==  1).OR.(nbondj == 2)) )  ztab(nlci-1,nlcj-1,:) = ptab(nlci-1,nlcj-1,:) ! East north
-            IF( (western_side) .AND. ((nbondj == -1).OR.(nbondj == 2)) )  ztab(2,2,:)           = ptab(2,2,:)           ! West south
-            IF( (western_side) .AND. ((nbondj ==  1).OR.(nbondj == 2)) )  ztab(2,nlcj-1,:)      = ptab(2,nlcj-1,:)      ! West north
-            
-            ! retrieve ice tracers
-            jm = 1
-            DO jl = 1, jpl
-               !
-               DO jj = j1, j2
-                  DO ji = i1, i2
-                     a_i (ji,jj,jl) = ztab(ji,jj,jm  ) * tmask(ji,jj,1)
-                     v_i (ji,jj,jl) = ztab(ji,jj,jm+1) * tmask(ji,jj,1)
-                     v_s (ji,jj,jl) = ztab(ji,jj,jm+2) * tmask(ji,jj,1)
-                     sv_i(ji,jj,jl) = ztab(ji,jj,jm+3) * tmask(ji,jj,1)
-                     oa_i(ji,jj,jl) = ztab(ji,jj,jm+4) * tmask(ji,jj,1)
-                     a_ip(ji,jj,jl) = ztab(ji,jj,jm+5) * tmask(ji,jj,1)
-                     v_ip(ji,jj,jl) = ztab(ji,jj,jm+6) * tmask(ji,jj,1)
-                     t_su(ji,jj,jl) = ztab(ji,jj,jm+7) * tmask(ji,jj,1)
-                  END DO
-               END DO
-               jm = jm + 8
-               !
-               DO jk = 1, nlay_s
-                  e_s(i1:i2,j1:j2,jk,jl) = ztab(i1:i2,j1:j2,jm) * tmask(i1:i2,j1:j2,1)
-                  jm = jm + 1
-               END DO
-               !
-               DO jk = 1, nlay_i
-                  e_i(i1:i2,j1:j2,jk,jl) = ztab(i1:i2,j1:j2,jm) * tmask(i1:i2,j1:j2,1)
-                  jm = jm + 1
-               END DO
-               !
-            END DO
-          
-         ENDIF  ! nbghostcells=1
+!!==> clem: this interpolation does not work because it creates negative values, due
+!           to negative coefficients when mixing points (for ex. z7)
+!!
+!         ELSE                          ! ==> complex interpolation (only one ghost cell available)
+!            !! Use a more complex interpolation since we mix solutions over a couple of grid points
+!            !! it is advised to use it for fields modified by high order schemes (e.g. advection UM5...)
+!            !! clem: for some reason (I don't know why), the following lines do not work 
+!            !        I think there is an issue with Agrif_SpecialValue here (not taken into account properly)
+!            ! record ztab
+!            jm = 1
+!            DO jl = 1, jpl
+!               ztab(:,:,jm  ) = a_i (:,:,jl)
+!               ztab(:,:,jm+1) = v_i (:,:,jl)
+!               ztab(:,:,jm+2) = v_s (:,:,jl)
+!               ztab(:,:,jm+3) = sv_i(:,:,jl)
+!               ztab(:,:,jm+4) = oa_i(:,:,jl)
+!               ztab(:,:,jm+5) = a_ip(:,:,jl)
+!               ztab(:,:,jm+6) = v_ip(:,:,jl)
+!               ztab(:,:,jm+7) = t_su(:,:,jl)
+!               jm = jm + 8
+!               DO jk = 1, nlay_s
+!                  ztab(:,:,jm) = e_s(:,:,jk,jl)
+!                  jm = jm + 1
+!               END DO
+!               DO jk = 1, nlay_i
+!                  ztab(:,:,jm) = e_i(:,:,jk,jl)
+!                  jm = jm + 1
+!               END DO
+!               !
+!            END DO
+!            !
+!            ! borders of the domain
+!            western_side  = (nb == 1).AND.(ndir == 1)  ;  eastern_side  = (nb == 1).AND.(ndir == 2)
+!            southern_side = (nb == 2).AND.(ndir == 1)  ;  northern_side = (nb == 2).AND.(ndir == 2)
+!            !
+!            ! spatial smoothing
+!            zrhox = Agrif_Rhox()
+!            z1 =      ( zrhox - 1. ) * 0.5 
+!            z3 =      ( zrhox - 1. ) / ( zrhox + 1. )
+!            z6 = 2. * ( zrhox - 1. ) / ( zrhox + 1. )
+!            z7 =    - ( zrhox - 1. ) / ( zrhox + 3. )
+!            z2 = 1. - z1
+!            z4 = 1. - z3
+!            z5 = 1. - z6 - z7
+!            !
+!            ! Remove corners
+!            imin = i1  ;  imax = i2  ;  jmin = j1  ;  jmax = j2
+!            IF( (nbondj == -1) .OR. (nbondj == 2) )   jmin = 3
+!            IF( (nbondj == +1) .OR. (nbondj == 2) )   jmax = nlcj-2
+!            IF( (nbondi == -1) .OR. (nbondi == 2) )   imin = 3
+!            IF( (nbondi == +1) .OR. (nbondi == 2) )   imax = nlci-2
+!
+!            ! smoothed fields
+!            IF( eastern_side ) THEN
+!               ztab(nlci,j1:j2,:) = z1 * ptab(nlci,j1:j2,:) + z2 * ptab(nlci-1,j1:j2,:)
+!               DO jj = jmin, jmax
+!                  rswitch = 0.
+!                  IF( u_ice(nlci-2,jj) > 0._wp ) rswitch = 1.
+!                  ztab(nlci-1,jj,:) = ( 1. - umask(nlci-2,jj,1) ) * ztab(nlci,jj,:)  &
+!                     &                +      umask(nlci-2,jj,1)   *  &
+!                     &                ( ( 1. - rswitch ) * ( z4 * ztab(nlci,jj,:)   + z3 * ztab(nlci-2,jj,:) )  &
+!                     &                  +      rswitch   * ( z6 * ztab(nlci-2,jj,:) + z5 * ztab(nlci,jj,:) + z7 * ztab(nlci-3,jj,:) ) )
+!                  ztab(nlci-1,jj,:) = ztab(nlci-1,jj,:) * tmask(nlci-1,jj,1)
+!               END DO
+!            ENDIF
+!            ! 
+!            IF( northern_side ) THEN
+!               ztab(i1:i2,nlcj,:) = z1 * ptab(i1:i2,nlcj,:) + z2 * ptab(i1:i2,nlcj-1,:)
+!               DO ji = imin, imax
+!                  rswitch = 0.
+!                  IF( v_ice(ji,nlcj-2) > 0._wp ) rswitch = 1.
+!                  ztab(ji,nlcj-1,:) = ( 1. - vmask(ji,nlcj-2,1) ) * ztab(ji,nlcj,:)  &
+!                     &                +      vmask(ji,nlcj-2,1)   *  &
+!                     &                ( ( 1. - rswitch ) * ( z4 * ztab(ji,nlcj,:)   + z3 * ztab(ji,nlcj-2,:) ) &
+!                     &                  +      rswitch   * ( z6 * ztab(ji,nlcj-2,:) + z5 * ztab(ji,nlcj,:) + z7 * ztab(ji,nlcj-3,:) ) )
+!                  ztab(ji,nlcj-1,:) = ztab(ji,nlcj-1,:) * tmask(ji,nlcj-1,1)
+!               END DO
+!            END IF
+!            !
+!            IF( western_side) THEN
+!               ztab(1,j1:j2,:) = z1 * ptab(1,j1:j2,:) + z2 * ptab(2,j1:j2,:)
+!               DO jj = jmin, jmax
+!                  rswitch = 0.
+!                  IF( u_ice(2,jj) < 0._wp ) rswitch = 1.
+!                  ztab(2,jj,:) = ( 1. - umask(2,jj,1) ) * ztab(1,jj,:)  &
+!                     &           +      umask(2,jj,1)   *   &
+!                     &           ( ( 1. - rswitch ) * ( z4 * ztab(1,jj,:) + z3 * ztab(3,jj,:) ) &
+!                     &             +      rswitch   * ( z6 * ztab(3,jj,:) + z5 * ztab(1,jj,:) + z7 * ztab(4,jj,:) ) )
+!                  ztab(2,jj,:) = ztab(2,jj,:) * tmask(2,jj,1)
+!               END DO
+!            ENDIF
+!            !
+!            IF( southern_side ) THEN
+!               ztab(i1:i2,1,:) = z1 * ptab(i1:i2,1,:) + z2 * ptab(i1:i2,2,:)
+!               DO ji = imin, imax
+!                  rswitch = 0.
+!                  IF( v_ice(ji,2) < 0._wp ) rswitch = 1.
+!                  ztab(ji,2,:) = ( 1. - vmask(ji,2,1) ) * ztab(ji,1,:)  &
+!                     &           +      vmask(ji,2,1)   *  &
+!                     &           ( ( 1. - rswitch ) * ( z4 * ztab(ji,1,:) + z3 * ztab(ji,3,:) ) &
+!                     &             +      rswitch   * ( z6 * ztab(ji,3,:) + z5 * ztab(ji,1,:) + z7 * ztab(ji,4,:) ) )
+!                  ztab(ji,2,:) = ztab(ji,2,:) * tmask(ji,2,1)
+!               END DO
+!            END IF
+!            !
+!            ! Treatment of corners
+!            IF( (eastern_side) .AND. ((nbondj == -1).OR.(nbondj == 2)) )  ztab(nlci-1,2,:)      = ptab(nlci-1,2,:)      ! East south
+!            IF( (eastern_side) .AND. ((nbondj ==  1).OR.(nbondj == 2)) )  ztab(nlci-1,nlcj-1,:) = ptab(nlci-1,nlcj-1,:) ! East north
+!            IF( (western_side) .AND. ((nbondj == -1).OR.(nbondj == 2)) )  ztab(2,2,:)           = ptab(2,2,:)           ! West south
+!            IF( (western_side) .AND. ((nbondj ==  1).OR.(nbondj == 2)) )  ztab(2,nlcj-1,:)      = ptab(2,nlcj-1,:)      ! West north
+!            
+!            ! retrieve ice tracers
+!            jm = 1
+!            DO jl = 1, jpl
+!               !
+!               DO jj = j1, j2
+!                  DO ji = i1, i2
+!                     a_i (ji,jj,jl) = ztab(ji,jj,jm  ) * tmask(ji,jj,1)
+!                     v_i (ji,jj,jl) = ztab(ji,jj,jm+1) * tmask(ji,jj,1)
+!                     v_s (ji,jj,jl) = ztab(ji,jj,jm+2) * tmask(ji,jj,1)
+!                     sv_i(ji,jj,jl) = ztab(ji,jj,jm+3) * tmask(ji,jj,1)
+!                     oa_i(ji,jj,jl) = ztab(ji,jj,jm+4) * tmask(ji,jj,1)
+!                     a_ip(ji,jj,jl) = ztab(ji,jj,jm+5) * tmask(ji,jj,1)
+!                     v_ip(ji,jj,jl) = ztab(ji,jj,jm+6) * tmask(ji,jj,1)
+!                     t_su(ji,jj,jl) = ztab(ji,jj,jm+7) * tmask(ji,jj,1)
+!                  END DO
+!               END DO
+!               jm = jm + 8
+!               !
+!               DO jk = 1, nlay_s
+!                  e_s(i1:i2,j1:j2,jk,jl) = ztab(i1:i2,j1:j2,jm) * tmask(i1:i2,j1:j2,1)
+!                  jm = jm + 1
+!               END DO
+!               !
+!               DO jk = 1, nlay_i
+!                  e_i(i1:i2,j1:j2,jk,jl) = ztab(i1:i2,j1:j2,jm) * tmask(i1:i2,j1:j2,1)
+!                  jm = jm + 1
+!               END DO
+!               !
+!            END DO
+!          
+!         ENDIF  ! nbghostcells=1
          
          ! integrated values
@@ -361,4 +364,13 @@
          et_s(i1:i2,j1:j2)  = SUM( SUM( e_s(i1:i2,j1:j2,:,:), dim=4 ), dim=3 )
          et_i(i1:i2,j1:j2)  = SUM( SUM( e_i(i1:i2,j1:j2,:,:), dim=4 ), dim=3 )
+
+         at_ip(i1:i2,j1:j2) = SUM( a_ip(i1:i2,j1:j2,:), dim=3 ) ! melt ponds
+         vt_ip(i1:i2,j1:j2) = SUM( v_ip(i1:i2,j1:j2,:), dim=3 )
+         !
+         ato_i(i1:i2,j1:j2) = 1._wp - at_i(i1:i2,j1:j2)         ! open water fraction  
+
+         DO jl = 1, jpl
+            WHERE( tmask(i1:i2,j1:j2,1) == 0._wp )   t_su(i1:i2,j1:j2,jl) = rt0   ! to avoid a division by 0 in sbcblk.F90
+         END DO
          !
       ENDIF
