Index: /branches/2012/dev_r3387_LOCEAN6_AGRIF_LIM/NEMOGCM/NEMO/NST_SRC/agrif_lim2_interp.F90
===================================================================
--- /branches/2012/dev_r3387_LOCEAN6_AGRIF_LIM/NEMOGCM/NEMO/NST_SRC/agrif_lim2_interp.F90	(revision 3479)
+++ /branches/2012/dev_r3387_LOCEAN6_AGRIF_LIM/NEMOGCM/NEMO/NST_SRC/agrif_lim2_interp.F90	(revision 3479)
@@ -0,0 +1,610 @@
+MODULE agrif_lim2_interp
+   !!======================================================================
+   !!                       ***  MODULE agrif_lim2_update ***
+   !! Nesting module :  update surface ocean boundary condition over ice
+   !!                   from a child grif
+   !! Sea-Ice model  :  LIM 2.0 Sea ice model time-stepping
+   !!======================================================================
+   !! History :  2.0   !  04-2008  (F. Dupont)  initial version
+   !!            3.4   !  09-2012  (R. Benshila, C. Herbaut) update and EVP
+   !!----------------------------------------------------------------------
+#if defined key_agrif && defined key_lim2
+   !!----------------------------------------------------------------------
+   !!   'key_lim2'  :                                 LIM 2.0 sea-ice model
+   !!   'key_agrif' :                                 AGRIF library
+   !!----------------------------------------------------------------------
+   !!   agrif_interp_lim2   : update sea-ice model on boundaries or total
+   !!                         sea-ice area
+   !!  agrif_rhg_lim2_load  : interpolcation of ice velocities using Agrif
+   !!  agrif_rhg_lim2       : sub-interpolation of ice velocities for both 
+   !!                         splitting time and sea-ice time step
+   !!  agrif_interp_u_ice   : atomic routine to interpolate u_ice 
+   !!  agrif_interp_u_ice   : atomic routine to interpolate v_ice 
+   !!  agrif_trp_lim2_load  : interpolcation of ice properties using Agrif
+   !!  agrif_trp_lim2       : sub-interpolation of ice properties for  
+   !!                         sea-ice time step
+   !!  agrif_interp_u_ice   : atomic routine to interpolate ice properties 
+   !!----------------------------------------------------------------------
+   USE par_oce
+   USE dom_oce
+   USE sbc_oce
+   USE ice_2
+   USE dom_ice_2
+   USE agrif_ice
+
+   IMPLICIT NONE
+   PRIVATE
+
+   PUBLIC agrif_rhg_lim2_load, agrif_rhg_lim2
+   PUBLIC agrif_trp_lim2_load, agrif_trp_lim2
+   PUBLIC interp_u_ice, interp_v_ice
+   PUBLIC interp_adv_ice
+
+   !!----------------------------------------------------------------------
+   !! NEMO/NST 3.4 , NEMO Consortium (2012)
+   !! $Id$
+   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
+   !!----------------------------------------------------------------------
+
+CONTAINS
+
+# if defined key_lim2_vp
+   SUBROUTINE agrif_rhg_lim2_load
+      !!-----------------------------------------------------------------------
+      !!              *** ROUTINE agrif_rhg_lim2_load ***
+      !!
+      !!  ** Method  : need a special routine for dealing with exchanging data
+      !! between the child and parent grid during ice step
+      !!
+      !!-----------------------------------------------------------------------
+      !
+      IF (Agrif_Root()) RETURN
+
+      Agrif_SpecialValue=0.
+      Agrif_UseSpecialValue = .FALSE.
+      u_ice_nst(:,:) = 0.
+      v_ice_nst(:,:) = 0.
+      CALL Agrif_Bc_variable( u_ice_nst, u_ice_id ,procname=interp_u_ice, calledweight=1. )
+      CALL Agrif_Bc_variable( v_ice_nst, v_ice_id ,procname=interp_v_ice, calledweight=1. )
+      Agrif_SpecialValue=0.
+      Agrif_UseSpecialValue = .FALSE.
+      !
+   END SUBROUTINE agrif_rhg_lim2_load
+
+
+   SUBROUTINE agrif_rhg_lim2(pu_n,pv_n)
+      !!-----------------------------------------------------------------------
+      !!                 *** ROUTINE agrif_rhg_lim2 ***
+      !!
+      !!  ** Method  : we feel the boundaries with values stored above
+      !!-----------------------------------------------------------------------
+      REAL(wp), DIMENSION(jpi,0:jpj+1), INTENT(inout) :: pu_n, pv_n
+      !!
+      REAL(wp) :: zrhox, zrhoy
+      INTEGER :: ji,jj
+      !!-----------------------------------------------------------------------
+      !
+      IF (Agrif_Root()) RETURN
+
+      zrhox = Agrif_Rhox()
+      zrhoy = Agrif_Rhoy()
+
+      IF((nbondi == -1).OR.(nbondi == 2)) THEN
+         DO jj=2,jpj
+            pu_n(3,jj) = u_ice_nst(3,jj)/(zrhoy*e2f(2,jj-1))*tmu(3,jj)
+         END DO
+         DO jj=2,jpj
+            pv_n(3,jj) = v_ice_nst(3,jj)/(zrhox*e1f(2,jj-1))*tmu(3,jj)
+         END DO
+      ENDIF
+
+      IF((nbondi == 1).OR.(nbondi == 2)) THEN
+         DO jj=2,jpj
+            pu_n(nlci-1,jj) = u_ice_nst(nlci-1,jj)/(zrhoy*e2f(nlci-2,jj-1))*tmu(nlci-1,jj)
+         END DO
+         DO jj=2,jpj
+            pv_n(nlci-1,jj) = v_ice_nst(nlci-1,jj)/(zrhox*e1f(nlci-2,jj-1))*tmu(nlci-1,jj)
+         END DO
+      ENDIF
+
+      IF((nbondj == -1).OR.(nbondj == 2)) THEN
+         DO ji=2,jpi
+            pv_n(ji,3) = v_ice_nst(ji,3)/(zrhox*e1f(ji-1,2))*tmu(ji,3)
+         END DO
+         DO ji=2,jpi
+            pu_n(ji,3) = u_ice_nst(ji,3)/(zrhoy*e2f(ji-1,2))*tmu(ji,3)
+         END DO
+      ENDIF
+
+      IF((nbondj == 1).OR.(nbondj == 2)) THEN
+         DO ji=2,jpi
+            pv_n(ji,nlcj-1) = v_ice_nst(ji,nlcj-1)/(zrhox*e1f(ji-1,nlcj-2))*tmu(ji,nlcj-1)
+         END DO
+         DO ji=2,jpi
+            pu_n(ji,nlcj-1) = u_ice_nst(ji,nlcj-1)/(zrhoy*e2f(ji-1,nlcj-2))*tmu(ji,nlcj-1)
+         END DO
+      ENDIF
+      !
+   END SUBROUTINE agrif_rhg_lim2
+
+#else
+   SUBROUTINE agrif_rhg_lim2_load
+      !!-----------------------------------------------------------------------
+      !!              *** ROUTINE agrif_rhg_lim2_load ***
+      !!
+      !!  ** Method  : need a special routine for dealing with exchanging data
+      !!  between the child and parent grid during ice step
+      !!               we interpolate and store the boundary if needed, ie if
+      !!  we are in inside a new parent ice time step
+      !!-----------------------------------------------------------------------
+      REAL(wp), DIMENSION(jpi,jpj) :: zuice, zvice
+      INTEGER :: ji,jj
+      REAL(wp) :: zrhox, zrhoy
+      !!-----------------------------------------------------------------------
+      !
+      IF (Agrif_Root()) RETURN
+
+      IF( lim_nbstep == 1. ) THEN
+         !
+         ! switch old values by hand
+         u_ice_oe(:,:,1) =  u_ice_oe(:,:,2)
+         v_ice_oe(:,:,1) =  v_ice_oe(:,:,2)
+         u_ice_sn(:,:,1) =  u_ice_sn(:,:,2)
+         v_ice_sn(:,:,1) =  v_ice_sn(:,:,2)
+         ! interpolation of boundaries (called weight prevents AGRIF interpolation)
+         Agrif_SpecialValue=-9999.
+         Agrif_UseSpecialValue = .TRUE.
+         zuice = 0.
+         zvice = 0.
+         CALL Agrif_Bc_variable(zuice,u_ice_id,procname=interp_u_ice, calledweight=1.)
+         CALL Agrif_Bc_variable(zvice,v_ice_id,procname=interp_v_ice, calledweight=1.)
+         Agrif_SpecialValue=0.
+         Agrif_UseSpecialValue = .FALSE.
+         !  
+         zrhox = agrif_rhox() ;    zrhoy = agrif_rhoy()      
+         zuice(:,:) =  zuice(:,:)/(zrhoy*e2u(:,:))*umask(:,:,1)
+         zvice(:,:) =  zvice(:,:)/(zrhox*e1v(:,:))*vmask(:,:,1)
+         ! fill  boundaries
+         DO jj = 1, jpj
+            DO ji = 1, 2
+               u_ice_oe(ji,  jj,2) = zuice(ji       ,jj) 
+               u_ice_oe(ji+2,jj,2) = zuice(nlci+ji-3,jj)
+            END DO
+         END DO
+         DO jj = 1, jpj
+            v_ice_oe(2,jj,2) = zvice(2     ,jj) 
+            v_ice_oe(4,jj,2) = zvice(nlci-1,jj)
+         END DO
+         DO ji = 1, jpi
+            u_ice_sn(ji,2,2) = zuice(ji,2     ) 
+            u_ice_sn(ji,4,2) = zuice(ji,nlcj-1)
+         END DO
+         DO jj = 1, 2
+            DO ji = 1, jpi
+               v_ice_sn(ji,jj  ,2) = zvice(ji,jj       ) 
+               v_ice_sn(ji,jj+2,2) = zvice(ji,nlcj+jj-3)
+            END DO
+         END DO
+         !
+      ENDIF
+      !
+   END SUBROUTINE agrif_rhg_lim2_load
+
+
+   SUBROUTINE agrif_rhg_lim2( kiter, kitermax, cd_type)
+      !!-----------------------------------------------------------------------
+      !!                 *** ROUTINE agrif_rhg_lim2  ***
+      !!
+      !!  ** Method  : simple call to atomic routines using stored values to
+      !!  fill the boundaries depending of the position of the point and 
+      !!  computing factor for time interpolation 
+      !!-----------------------------------------------------------------------
+      INTEGER, INTENT(in) :: kiter, kitermax
+      CHARACTER(len=1), INTENT( in ) :: cd_type
+      !!    
+      REAL(wp) :: zalpha, zbeta
+      !!-----------------------------------------------------------------------
+      !
+      IF (Agrif_Root()) RETURN
+      zalpha = REAL(lim_nbstep,wp) / (Agrif_Rhot()*Agrif_PArent(nn_fsbc)/REAL(nn_fsbc))
+      zbeta  = REAL(kiter,wp) / kitermax
+      zbeta = zalpha * zbeta
+      SELECT CASE(cd_type)
+      CASE('U')
+         CALL ParcoursU( zbeta )
+      CASE('V')
+         CALL ParcoursV( zbeta )
+      END SELECT
+      !
+   END SUBROUTINE agrif_rhg_lim2
+
+
+   SUBROUTINE ParcoursU( pbeta )
+      !!-----------------------------------------------------------------------
+      !!                    *** ROUTINE parcoursU ***
+      !!
+      !!  ** Method  : time and spatial interpolation for U-point using values
+      !!  interpolated from the coarse grid and inside dvalues     
+      !!-----------------------------------------------------------------------
+      REAL(wp), INTENT(in) :: pbeta
+      !!
+      INTEGER :: ji, jj
+      !!-----------------------------------------------------------------------
+      !
+      IF((nbondi == -1).OR.(nbondi == 2)) THEN
+         DO jj=1,jpj
+            DO ji=1,2
+               u_ice(ji,jj) = (1-pbeta) * u_ice_oe(ji,jj,1) + pbeta * u_ice_oe(ji,jj,2)
+            END DO
+         END DO
+         DO jj=1,jpj
+            u_ice(2,jj) = 0.25*(u_ice(1,jj)+2.*u_ice(2,jj)+u_ice(3,jj))
+            u_ice(2,jj) = u_ice(2,jj) * umask(2,jj,1)
+         END DO
+      ENDIF
+
+      IF((nbondi == 1).OR.(nbondi == 2)) THEN
+         DO jj=1,jpj
+            DO ji=1,2
+               u_ice(nlci+ji-3,jj) = (1-pbeta) * u_ice_oe(ji+2,jj,1) + pbeta * u_ice_oe(ji+2,jj,2)
+            END DO
+         END DO
+         DO jj=1,jpj
+            u_ice(nlci-2,jj) = 0.25*(u_ice(nlci-3,jj)+2.*u_ice(nlci-2,jj)+u_ice(nlci-1,jj))
+            u_ice(nlci-2,jj) = u_ice(nlci-2,jj) * umask(nlci-2,jj,1)
+         END DO
+      ENDIF
+
+      IF((nbondj == -1).OR.(nbondj == 2)) THEN
+         DO ji=1,jpi
+            u_ice(ji,2) = (1-pbeta) * u_ice_sn(ji,2,1) + pbeta * u_ice_sn(ji,2,2)
+            u_ice(ji,2) = u_ice(ji,2)*umask(ji,2,1)
+         END DO
+      ENDIF
+
+      IF((nbondj == 1).OR.(nbondj == 2)) THEN
+         DO ji=1,jpi
+            u_ice(ji,nlcj-1) = (1-pbeta) * u_ice_sn(ji,4,1) + pbeta * u_ice_sn(ji,4,2)
+            u_ice(ji,nlcj-1) = u_ice(ji,nlcj-1)*umask(ji,nlcj-1,1)
+         END DO
+      ENDIF
+      !
+   END SUBROUTINE ParcoursU
+
+
+   SUBROUTINE ParcoursV( pbeta )
+      !!-----------------------------------------------------------------------
+      !!                    *** ROUTINE parcoursV ***
+      !! 
+      !!  ** Method  : time and spatial interpolation for V-point using values
+      !!  interpolated from the coarse grid and inside dvalues     
+      !!-----------------------------------------------------------------------
+      REAL(wp), INTENT(in) :: pbeta
+      !!
+      INTEGER :: ji, jj
+      !!-----------------------------------------------------------------------
+      !
+      IF((nbondi == -1).OR.(nbondi == 2)) THEN
+         DO jj=1,jpj
+            v_ice(2,jj) = (1-pbeta) * v_ice_oe(2,jj,1) + pbeta * v_ice_oe(2,jj,2)
+            v_ice(2,jj) = v_ice(2,jj) * vmask(2,jj,1)
+         END DO
+      ENDIF
+
+      IF((nbondi == 1).OR.(nbondi == 2)) THEN
+         DO jj=1,jpj
+            v_ice(nlci-1,jj) = (1-pbeta) * v_ice_oe(4,jj,1) + pbeta * v_ice_oe(4,jj,2)
+            v_ice(nlci-1,jj) = v_ice(nlci-1,jj)*vmask(nlci-1,jj,1)
+         END DO
+      ENDIF
+
+      IF((nbondj == -1).OR.(nbondj == 2)) THEN
+         DO jj=1,2
+            DO ji=1,jpi
+               v_ice(ji,jj) = (1-pbeta) * v_ice_sn(ji,jj,1) + pbeta * v_ice_sn(ji,jj,2)
+            END DO
+         END DO
+         DO ji=1,jpi
+            v_ice(ji,2)=0.25*(v_ice(ji,1)+2.*v_ice(ji,2)+v_ice(ji,3))
+            v_ice(ji,2)=v_ice(ji,2)*vmask(ji,2,1)
+         END DO
+      ENDIF
+
+      IF((nbondj == 1).OR.(nbondj == 2)) THEN
+         DO jj=1,2
+            DO ji=1,jpi
+               v_ice(ji,nlcj+jj-3) = (1-pbeta) * v_ice_sn(ji,jj+2,1) + pbeta * v_ice_sn(ji,jj+2,2)
+            END DO
+         END DO
+         DO ji=1,jpi
+            v_ice(ji,nlcj-2)=0.25*(v_ice(ji,nlcj-3)+2.*v_ice(ji,nlcj-2)+v_ice(ji,nlcj-1))
+            v_ice(ji,nlcj-2) = v_ice(ji,nlcj-2) * vmask(ji,nlcj-2,1)
+         END DO
+      ENDIF
+      !
+   END SUBROUTINE ParcoursV
+# endif
+   SUBROUTINE agrif_trp_lim2_load
+      !!-----------------------------------------------------------------------
+      !!                 *** ROUTINE agrif_trp_lim2_load ***
+      !!
+      !!  ** Method  : need a special routine for dealing with exchanging data
+      !!  between the child and parent grid during ice step
+      !!               we interpolate and store the boundary if needed, ie if
+      !!  we are in inside a new parent ice time step
+     !!-----------------------------------------------------------------------
+      REAL(wp), DIMENSION(jpi,jpj,7) :: ztab 
+      INTEGER :: ji,jj,jn
+      !!-----------------------------------------------------------------------
+      !
+      IF (Agrif_Root()) RETURN
+      IF( lim_nbstep == 1. ) THEN
+         !
+         ! switch old values
+         adv_ice_oe(:,:,:,1) =  adv_ice_oe(:,:,:,2)
+         adv_ice_sn(:,:,:,1) =  adv_ice_sn(:,:,:,2)
+         ! interpolation of boundaries
+         ztab(:,:,:) = 0.
+         Agrif_SpecialValue=-9999.
+         Agrif_UseSpecialValue = .TRUE.
+         CALL Agrif_Bc_variable( ztab, adv_ice_id ,procname=interp_adv_ice,calledweight=1. )
+         Agrif_SpecialValue=0.
+         Agrif_UseSpecialValue = .FALSE.
+         !  
+         ! fill  boundaries
+         DO jn =1,7
+            DO jj = 1, jpj
+               DO ji=1,2
+                  adv_ice_oe(ji  ,jj,jn,2) = ztab(ji       ,jj,jn) 
+                  adv_ice_oe(ji+2,jj,jn,2) = ztab(nlci-2+ji,jj,jn)
+               END DO
+            END DO
+         END DO
+
+         Do jn =1,7
+            Do jj =1,2
+               DO ji = 1, jpi
+                  adv_ice_sn(ji,jj  ,jn,2) = ztab(ji,jj       ,jn) 
+                  adv_ice_sn(ji,jj+2,jn,2) = ztab(ji,nlcj-2+jj,jn)
+               END DO
+            END DO
+         END DO
+         !
+      ENDIF
+      !
+   END SUBROUTINE agrif_trp_lim2_load
+
+
+   SUBROUTINE agrif_trp_lim2
+      !!-----------------------------------------------------------------------
+      !!                  *** ROUTINE agrif_trp_lim2 ***
+      !!
+      !!  ** Method  : time coefficient and call to atomic routines
+      !!-----------------------------------------------------------------------
+      INTEGER :: ji,jj,jn
+      REAL(wp) :: zalpha
+      REAL(wp), DIMENSION(jpi,jpj,7) :: ztab 
+      !!-----------------------------------------------------------------------      
+      !
+      IF (Agrif_Root()) RETURN
+
+      zalpha = REAL(lim_nbstep,wp) / (Agrif_Rhot()*Agrif_PArent(nn_fsbc)/REAL(nn_fsbc))
+      !
+      ztab(:,:,:) = 0.e0
+      DO jn =1,7
+         DO jj =1,2
+            DO ji = 1, jpi
+               ztab(ji,jj        ,jn) = (1-zalpha)*adv_ice_sn(ji,jj  ,jn,1) + zalpha*adv_ice_sn(ji,jj  ,jn,2) 
+               ztab(ji,nlcj-2+jj ,jn) = (1-zalpha)*adv_ice_sn(ji,jj+2,jn,1) + zalpha*adv_ice_sn(ji,jj+2,jn,2) 
+            END DO
+         END DO
+      END DO
+
+      DO jn =1,7
+         DO jj = 1, jpj
+            DO ji=1,2
+               ztab(ji       ,jj,jn) = (1-zalpha)*adv_ice_oe(ji  ,jj,jn,1) + zalpha*adv_ice_oe(ji  ,jj,jn,2) 
+               ztab(nlci-2+ji,jj,jn) = (1-zalpha)*adv_ice_oe(ji+2,jj,jn,1) + zalpha*adv_ice_oe(ji+2,jj,jn,2) 
+            END DO
+         END DO
+      END DO
+      !
+      CALL parcoursT( ztab(:,:, 1), frld  )
+      CALL parcoursT( ztab(:,:, 2), hicif )
+      CALL parcoursT( ztab(:,:, 3), hsnif )
+      CALL parcoursT( ztab(:,:, 4), tbif(:,:,1) )
+      CALL parcoursT( ztab(:,:, 5), tbif(:,:,2) )
+      CALL parcoursT( ztab(:,:, 6), tbif(:,:,3) )
+      CALL parcoursT( ztab(:,:, 7), qstoif )
+      !
+   END SUBROUTINE agrif_trp_lim2
+
+
+   SUBROUTINE parcoursT ( pinterp, pfinal )
+      !!-----------------------------------------------------------------------
+      !!                    *** ROUTINE parcoursT ***
+      !!
+      !!  ** Method  : fill boundaries for T points
+      !!-----------------------------------------------------------------------
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(in)    :: pinterp
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: pfinal
+      !!
+      REAL(wp) :: zbound, zvbord
+      REAL(wp), DIMENSION(jpi,jpj) ::  zui_u, zvi_v
+      INTEGER :: ji, jj
+      !!-----------------------------------------------------------------------
+      !
+      zui_u = 0.e0
+      zvi_v = 0.e0
+      ! zvbord factor between 1 and 2 to take into account slip or no-slip boundary conditions.
+      zbound=0.
+      zvbord = 1.0 + ( 1.0 - zbound )
+#if defined key_lim2_vp
+      DO jj = 1, jpjm1
+         DO ji = 1, jpim1
+            zui_u(ji,jj) = ( u_ice(ji+1,jj  ) + u_ice(ji+1,jj+1) ) / ( MAX( tmu(ji+1,jj  ) + tmu(ji+1,jj+1), zvbord ) )
+            zvi_v(ji,jj) = ( v_ice(ji  ,jj+1) + v_ice(ji+1,jj+1) ) / ( MAX( tmu(ji  ,jj+1) + tmu(ji+1,jj+1), zvbord ) )
+         END DO
+      END DO
+#else
+      zui_u(:,:) = u_ice(:,:)
+      zvi_v(:,:) = v_ice(:,:)
+#endif
+
+      IF((nbondi == -1).OR.(nbondi == 2)) THEN
+         DO jj=1,jpj
+            !            IF (zui_u(2,jj).EQ.0.) THEN
+            !               pfinal (2,jj) = pfinal (1,jj) * tms(2,jj)
+            !            ELSE
+            pfinal(2,jj) = 0.25* pinterp(1,jj) + 0.5 *  pinterp(2,jj) + 0.25 *pfinal(3,jj)
+            !            ENDIF
+         END DO
+      ENDIF
+
+      IF((nbondj == -1).OR.(nbondj == 2)) THEN
+         DO ji=1,jpi
+            !            IF (zvi_v(ji,2).EQ.0.) THEN
+            !               pfinal (ji,2) = pfinal (ji,1) * tms(ji,2)
+            !            ELSE
+            pfinal(ji,2) = 0.25* pinterp(ji,1) + 0.5 *  pinterp(ji,2) + 0.25 *pfinal(ji,3)
+            !            ENDIF 
+         END DO
+      ENDIF
+
+
+      IF((nbondi == 1).OR.(nbondi == 2)) THEN
+         DO jj=1,jpj
+            !            IF (zui_u(nlci-2,jj).EQ.0.) THEN
+            !               pfinal(nlci-1,jj) = pfinal (nlci,jj) * tms(nlci-1,jj)
+            !            ELSE
+            pfinal(nlci-1,jj) = 0.25* pinterp(nlci,jj) + 0.5 *  pinterp(nlci-1,jj) + 0.25 *pfinal(nlci-2,jj)
+            !           ENDIF
+         END DO
+      ENDIF
+
+      IF((nbondj == 1).OR.(nbondj == 2)) THEN
+         DO ji=1,jpi
+            !            IF (zvi_v(ji,nlcj-2).EQ.0.) THEN
+            !               pfinal (ji,nlcj-1) =  pfinal(ji,nlcj) * tms(ji,nlcj-1)
+            !            ELSE
+            pfinal(ji,nlcj-1) = 0.25* pinterp(ji,nlcj) + 0.5 *  pinterp(ji,nlcj-1) + 0.25 *pfinal(ji,nlcj-2)
+            !            ENDIF
+         END DO
+      ENDIF
+
+
+      pfinal (:,:) = pfinal (:,:) * tms(:,:)
+      ! 
+   END SUBROUTINE parcoursT
+
+
+   SUBROUTINE interp_u_ice( tabres, i1, i2, j1, j2 )
+      !!-----------------------------------------------------------------------
+      !!                     *** ROUTINE interp_u_ice ***
+      !!-----------------------------------------------------------------------
+      INTEGER, INTENT(in) :: i1, i2, j1, j2
+      REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) :: tabres
+      !!
+      INTEGER :: ji,jj
+      !!-----------------------------------------------------------------------
+      !
+#if defined key_lim2_vp
+      DO jj=MAX(j1,2),j2
+         DO ji=MAX(i1,2),i2
+            IF( tmu(ji,jj) == 0. ) THEN
+               tabres(ji,jj) = -9999.
+            ELSE
+               tabres(ji,jj) = e2f(ji-1,jj-1) * u_ice(ji,jj)
+            ENDIF
+         END DO
+      END DO
+#else
+      DO jj= j1, j2
+         DO ji= i1, i2
+            IF( umask(ji,jj,1) == 0. ) THEN
+               tabres(ji,jj) = -9999.
+            ELSE
+               tabres(ji,jj) = e2u(ji,jj) * u_ice(ji,jj)
+            ENDIF
+         END DO
+      END DO
+#endif
+   END SUBROUTINE interp_u_ice
+
+
+   SUBROUTINE interp_v_ice( tabres, i1, i2, j1, j2 )
+      !!-----------------------------------------------------------------------
+      !!                    *** ROUTINE interp_v_ice ***
+      !!-----------------------------------------------------------------------      
+      INTEGER, INTENT(in) :: i1, i2, j1, j2
+      REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) :: tabres
+      !!
+      INTEGER :: ji, jj
+      !!-----------------------------------------------------------------------
+      !
+#if defined key_lim2_vp
+      DO jj=MAX(j1,2),j2
+         DO ji=MAX(i1,2),i2
+            IF( tmu(ji,jj) == 0. ) THEN
+               tabres(ji,jj) = -9999.
+            ELSE
+               tabres(ji,jj) = e1f(ji-1,jj-1) * v_ice(ji,jj)
+            ENDIF
+         END DO
+      END DO
+#else
+      DO jj= j1 ,j2
+         DO ji = i1, i2
+            IF( vmask(ji,jj,1) == 0. ) THEN
+               tabres(ji,jj) = -9999.
+            ELSE
+               tabres(ji,jj) = e1v(ji,jj) * v_ice(ji,jj)
+            ENDIF
+         END DO
+      END DO
+#endif
+   END SUBROUTINE interp_v_ice
+
+
+   SUBROUTINE interp_adv_ice( tabres, i1, i2, j1, j2 )
+      !!-----------------------------------------------------------------------
+      !!                    *** ROUTINE interp_adv_ice ***                           
+      !!
+      !! ** Purpose : fill an array with  ice variables
+      !!              to be advected
+      !!              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 :: ji, jj, jk
+      !!-----------------------------------------------------------------------
+      !
+      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
+      !
+   END SUBROUTINE interp_adv_ice
+
+#else
+CONTAINS
+   SUBROUTINE agrif_lim2_interp_empty
+      !!---------------------------------------------
+      !!   *** ROUTINE agrif_lim2_interp_empty ***
+      !!---------------------------------------------
+      WRITE(*,*)  'agrif_lim2_interp : You should not have seen this print! error?'
+   END SUBROUTINE agrif_lim2_interp_empty
+#endif
+END MODULE agrif_lim2_interp
