Index: /trunk/AGRIF/AGRIF_FILES/modbc.F
===================================================================
--- /trunk/AGRIF/AGRIF_FILES/modbc.F	(revision 778)
+++ /trunk/AGRIF/AGRIF_FILES/modbc.F	(revision 779)
@@ -1046,6 +1046,6 @@
       if (child % var % interpIndex 
      &        /= Agrif_Curgrid % parent % ngridstep ) then
-       child%var%oldvalues2d(kindex:kindex+newsize-1,1)=
-     &        child%var%oldvalues2d(kindex:kindex+newsize-1,2)
+       child%var%oldvalues2d(1,kindex:kindex+newsize-1)=
+     &        child%var%oldvalues2d(2,kindex:kindex+newsize-1)
       endif
 
@@ -1053,6 +1053,7 @@
       CASE (1)
 
+!CDIR ALTCODE
          do ir=bounds(1,1),bounds(1,2)
-            child%var%oldvalues2d(kindex,2) =
+            child%var%oldvalues2d(2,kindex) =
      &           child%var%array1(ir)
             kindex = kindex + 1
@@ -1062,6 +1063,7 @@
 
         do jr=bounds(2,1),bounds(2,2)
+!CDIR ALTCODE
            do ir=bounds(1,1),bounds(1,2)
-            child%var%oldvalues2d(kindex,2) =
+            child%var%oldvalues2d(2,kindex) =
      &           child%var%array2(ir,jr)
             kindex = kindex + 1
@@ -1072,6 +1074,7 @@
         do kr=bounds(3,1),bounds(3,2)
           do jr=bounds(2,1),bounds(2,2)
+!CDIR ALTCODE
              do ir=bounds(1,1),bounds(1,2)
-            child%var%oldvalues2d(kindex,2) =
+            child%var%oldvalues2d(2,kindex) =
      &           child%var%array3(ir,jr,kr)
             kindex = kindex + 1
@@ -1084,6 +1087,7 @@
            do kr=bounds(3,1),bounds(3,2)
              do jr=bounds(2,1),bounds(2,2)
+!CDIR ALTCODE
                do ir=bounds(1,1),bounds(1,2)
-            child%var%oldvalues2d(kindex,2) =
+            child%var%oldvalues2d(2,kindex) =
      &           child%var%array4(ir,jr,kr,lr)
             kindex = kindex + 1
@@ -1098,6 +1102,7 @@
             do kr=bounds(3,1),bounds(3,2) 
               do jr=bounds(2,1),bounds(2,2)
+!CDIR ALTCODE
                  do ir=bounds(1,1),bounds(1,2)
-            child%var%oldvalues2d(kindex,2) =
+            child%var%oldvalues2d(2,kindex) =
      &           child%var%array5(ir,jr,kr,lr,mr)
             kindex = kindex + 1
@@ -1114,6 +1119,7 @@
                do kr=bounds(3,1),bounds(3,2)
                 do jr=bounds(2,1),bounds(2,2)
+!CDIR ALTCODE
                    do ir=bounds(1,1),bounds(1,2)
-            child%var%oldvalues2d(kindex,2) =
+            child%var%oldvalues2d(2,kindex) =
      &           child%var%array6(ir,jr,kr,lr,mr,nr)
             kindex = kindex + 1
@@ -1165,8 +1171,9 @@
       CASE (1)
 
+!CDIR ALTCODE
          do ir=bounds(1,1),bounds(1,2)
                 child%var%array1(ir) =
-     &           c2t*child % var % oldvalues2d(kindex,1)   
-     &         + c1t*child % var % oldvalues2d(kindex,2)    
+     &           c2t*child % var % oldvalues2d(1,kindex)   
+     &         + c1t*child % var % oldvalues2d(2,kindex)    
             kindex = kindex + 1
          enddo       
@@ -1175,8 +1182,9 @@
 
         do jr=bounds(2,1),bounds(2,2)
+!CDIR ALTCODE
            do ir=bounds(1,1),bounds(1,2)
                 child%var%array2(ir,jr) =
-     &           c2t*child % var % oldvalues2d(kindex,1)   
-     &         + c1t*child % var % oldvalues2d(kindex,2) 
+     &           c2t*child % var % oldvalues2d(1,kindex)   
+     &         + c1t*child % var % oldvalues2d(2,kindex) 
             kindex = kindex + 1
            enddo
@@ -1186,8 +1194,9 @@
         do kr=bounds(3,1),bounds(3,2)
           do jr=bounds(2,1),bounds(2,2)
+!CDIR ALTCODE
              do ir=bounds(1,1),bounds(1,2)
                 child%var%array3(ir,jr,kr) =
-     &           c2t*child % var % oldvalues2d(kindex,1)   
-     &         + c1t*child % var % oldvalues2d(kindex,2) 
+     &           c2t*child % var % oldvalues2d(1,kindex)   
+     &         + c1t*child % var % oldvalues2d(2,kindex) 
             kindex = kindex + 1
              enddo
@@ -1199,8 +1208,9 @@
            do kr=bounds(3,1),bounds(3,2)
              do jr=bounds(2,1),bounds(2,2)
+!CDIR ALTCODE
                do ir=bounds(1,1),bounds(1,2)
                 child%var%array4(ir,jr,kr,lr) =
-     &           c2t*child % var % oldvalues2d(kindex,1)   
-     &         + c1t*child % var % oldvalues2d(kindex,2) 
+     &           c2t*child % var % oldvalues2d(1,kindex)   
+     &         + c1t*child % var % oldvalues2d(2,kindex) 
             kindex = kindex + 1
                enddo
@@ -1214,8 +1224,9 @@
             do kr=bounds(3,1),bounds(3,2) 
               do jr=bounds(2,1),bounds(2,2)
+!CDIR ALTCODE
                  do ir=bounds(1,1),bounds(1,2)
                 child%var%array5(ir,jr,kr,lr,mr) =
-     &           c2t*child % var % oldvalues2d(kindex,1)   
-     &         + c1t*child % var % oldvalues2d(kindex,2) 
+     &           c2t*child % var % oldvalues2d(1,kindex)   
+     &         + c1t*child % var % oldvalues2d(2,kindex) 
             kindex = kindex + 1
                  enddo
@@ -1231,8 +1242,9 @@
                do kr=bounds(3,1),bounds(3,2)
                 do jr=bounds(2,1),bounds(2,2)
+!CDIR ALTCODE
                    do ir=bounds(1,1),bounds(1,2)
                 child%var%array6(ir,jr,kr,lr,mr,nr) =
-     &           c2t*child % var % oldvalues2d(kindex,1)   
-     &         + c1t*child % var % oldvalues2d(kindex,2) 
+     &           c2t*child % var % oldvalues2d(1,kindex)   
+     &         + c1t*child % var % oldvalues2d(2,kindex) 
             kindex = kindex + 1
                    enddo
@@ -1277,5 +1289,5 @@
       if (.NOT. associated(child % var % oldvalues2d)) then
 C
-          allocate(child % var % oldvalues2d(newsize,2))
+          allocate(child % var % oldvalues2d(2,newsize))
 C  
           child % var % oldvalues2d=0.
@@ -1283,8 +1295,8 @@
         else
 C
-          if (SIZE(child % var % oldvalues2d,1) < newsize) then   
-C
-              allocate(tempoldvalues(SIZE(child % var %
-     &                                    oldvalues2d,1),2))
+          if (SIZE(child % var % oldvalues2d,2) < newsize) then   
+C
+              allocate(tempoldvalues(2,SIZE(child % var %
+     &                                    oldvalues2d,2)))
 C
               tempoldvalues = child % var % oldvalues2d
@@ -1292,9 +1304,9 @@
               deallocate(child % var % oldvalues2d)
 C
-              allocate(child % var % oldvalues2d(newsize,2))
+              allocate(child % var % oldvalues2d(2,newsize))
 C            
               child%var%oldvalues2d=0.
 C
-              child % var % oldvalues2d(1:SIZE(tempoldvalues,1),:) = 
+              child % var % oldvalues2d(:,1:SIZE(tempoldvalues,1)) = 
      &        tempoldvalues(:,:)
 C
Index: /trunk/AGRIF/AGRIF_FILES/modbcfunction.F
===================================================================
--- /trunk/AGRIF/AGRIF_FILES/modbcfunction.F	(revision 778)
+++ /trunk/AGRIF/AGRIF_FILES/modbcfunction.F	(revision 779)
@@ -268,5 +268,5 @@
 
         Allocate(
-     &    Agrif_Curgrid%tabvars(tabvarsindic)%var % oldvalues2D(1,2))
+     &    Agrif_Curgrid%tabvars(tabvarsindic)%var % oldvalues2D(2,1))
           Agrif_Curgrid%tabvars(tabvarsindic)%var % oldvalues2D = 0. 
        ENDIF      
Index: /trunk/AGRIF/AGRIF_FILES/modinterp.F
===================================================================
--- /trunk/AGRIF/AGRIF_FILES/modinterp.F	(revision 778)
+++ /trunk/AGRIF/AGRIF_FILES/modinterp.F	(revision 779)
@@ -42,4 +42,5 @@
 C      
       IMPLICIT NONE
+      logical,  private:: precomputedone(7) = .FALSE.
 C
       CONTAINS
@@ -1351,4 +1352,5 @@
      &                  s_parent(nbdim),s_child(nbdim),
      &                  ds_parent(nbdim),ds_child(nbdim),coeffraf)
+                
 C                
       Return
@@ -1399,4 +1401,11 @@
       INTEGER i,j
       INTEGER :: coeffraf
+      REAL   , DIMENSION(
+     &                pttab_child(nbdim):petab_child(nbdim),
+     &                pttab_child(nbdim-1):petab_child(nbdim-1)
+     &                ) :: tabout_trsp
+      REAL, DIMENSION(indmin(nbdim):indmax(nbdim),
+     &        pttab_child(nbdim-1):petab_child(nbdim-1)) :: tabtemp_trsp
+
 C
 C
@@ -1405,6 +1414,35 @@
 C     Commentaire perso : nbdim vaut toujours 2 ici.
 C
+      coeffraf = nint ( ds_parent(1) / ds_child(1) )
+      IF((TypeInterp(1) == Agrif_Linear) .AND. (coeffraf /= 1 ) )THEN
+
+!---CDIR NEXPAND
+          IF(.NOT. precomputedone(1) ) call linear1Dprecompute2D(
+     &          indmax(2)-indmin(2)+1,   
+     &          indmax(1)-indmin(1)+1,   
+     &          petab_child(1)-pttab_child(1)+1,
+     &          s_parent(1),s_child(1),ds_parent(1),ds_child(1),1)
+!---CDIR NEXPAND
+           call linear1daftercompute(tabin,tabtemp,
+     &          size(tabin), size(tabtemp),  
+     &          s_parent(1),s_child(1),ds_parent(1),ds_child(1),1)
+ 
+      ELSEIF((TypeInterp(1) == Agrif_PPM) .AND. (coeffraf /= 1 ) )THEN
+!---CDIR NEXPAND
+          IF(.NOT. precomputedone(1) ) call ppm1Dprecompute2D(
+     &          indmax(2)-indmin(2)+1,   
+     &          indmax(1)-indmin(1)+1,   
+     &          petab_child(1)-pttab_child(1)+1,
+     &          s_parent(1),s_child(1),ds_parent(1),ds_child(1),1)
+!---CDIR NEXPAND
+           call ppm1daftercompute(tabin,tabtemp,
+     &          size(tabin), size(tabtemp),  
+     &          s_parent(1),s_child(1),ds_parent(1),ds_child(1),1)
+
+      ELSE
+
       do j = indmin(nbdim),indmax(nbdim)
 C        
+!---CDIR NEXPAND
         Call Agrif_Interp_1D_recursive(TypeInterp(1),
      &         tabin(indmin(nbdim-1):indmax(nbdim-1),j),
@@ -1416,18 +1454,53 @@
 C        
       enddo
-C        
+      ENDIF   
+
       coeffraf = nint(ds_parent(nbdim)/ds_child(nbdim))
       
+      tabtemp_trsp = TRANSPOSE(tabtemp)
+
+      IF((TypeInterp(2) == Agrif_Linear) .AND. (coeffraf /= 1 ) )THEN
+
+!---CDIR NEXPAND
+          IF(.NOT. precomputedone(2) ) call linear1Dprecompute2D(
+     &          petab_child(1)-pttab_child(1)+1,
+     &          indmax(2)-indmin(2)+1,
+     &          petab_child(2)-pttab_child(2)+1,
+     &          s_parent(2),s_child(2),ds_parent(2),ds_child(2),2)
+!---CDIR NEXPAND
+           call linear1daftercompute(tabtemp_trsp,tabout_trsp,
+     &          size(tabtemp_trsp), size(tabout_trsp),
+     &          s_parent(2),s_child(2),ds_parent(2),ds_child(2),2)
+
+      ELSEIF((TypeInterp(2) == Agrif_PPM) .AND. (coeffraf /= 1 ) )THEN
+
+!---CDIR NEXPAND
+           IF(.NOT. precomputedone(1) )call ppm1Dprecompute2D(
+     &          petab_child(1)-pttab_child(1)+1,
+     &          indmax(2)-indmin(2)+1,
+     &          petab_child(2)-pttab_child(2)+1,
+     &          s_parent(2),s_child(2),ds_parent(2),ds_child(2),2)
+!---CDIR NEXPAND
+           call ppm1daftercompute(tabtemp_trsp,tabout_trsp,
+     &          size(tabtemp_trsp), size(tabout_trsp),
+     &          s_parent(2),s_child(2),ds_parent(2),ds_child(2),2)
+
+      ELSE
       do i=pttab_child(nbdim-1),petab_child(nbdim-1)
 C
+!---CDIR NEXPAND
         Call Agrif_InterpBase(TypeInterp(2),
-     &           tabtemp(i,indmin(nbdim):indmax(nbdim)),
-     &                  tabout(i,pttab_child(nbdim):petab_child(nbdim)),
+     &           tabtemp_trsp(indmin(nbdim):indmax(nbdim),i),
+     &           tabout_trsp(pttab_child(nbdim):petab_child(nbdim),i),
      &           indmin(nbdim),indmax(nbdim),
      &           pttab_child(nbdim),petab_child(nbdim),
      &           s_parent(nbdim),s_child(nbdim),
      &           ds_parent(nbdim),ds_child(nbdim),coeffraf)
+
 C        
       enddo
+      ENDIF
+      
+      tabout = TRANSPOSE(tabout_trsp)
 C
       Return
@@ -1476,4 +1549,38 @@
 C
 C
+      coeffraf = nint ( ds_parent(1) / ds_child(1) )
+      IF((TypeInterp(1) == Agrif_Linear) .AND. (coeffraf/=1))THEN 
+        Call linear1Dprecompute2D(
+     &          indmax(2)-indmin(2)+1,
+     &          indmax(1)-indmin(1)+1,
+     &          petab_child(1)-pttab_child(1)+1,
+     &          s_parent(1),s_child(1),ds_parent(1),ds_child(1),1)
+      precomputedone(1) = .TRUE. 
+      ELSEIF((TypeInterp(1) == Agrif_PPM) .AND. (coeffraf/=1))THEN
+        Call ppm1Dprecompute2D(
+     &          indmax(2)-indmin(2)+1,
+     &          indmax(1)-indmin(1)+1,
+     &          petab_child(1)-pttab_child(1)+1,
+     &          s_parent(1),s_child(1),ds_parent(1),ds_child(1),1)
+      precomputedone(1) = .TRUE.
+      ENDIF
+
+      coeffraf = nint ( ds_parent(2) / ds_child(2) )
+      IF((TypeInterp(2) == Agrif_Linear) .AND. (coeffraf/=1)) THEN
+         Call linear1Dprecompute2D(
+     &          petab_child(1)-pttab_child(1)+1,
+     &          indmax(2)-indmin(2)+1,
+     &          petab_child(2)-pttab_child(2)+1,
+     &          s_parent(2),s_child(2),ds_parent(2),ds_child(2),2)
+      precomputedone(2) = .TRUE. 
+      ELSEIF((TypeInterp(2) == Agrif_PPM) .AND. (coeffraf/=1)) THEN
+         Call ppm1Dprecompute2D(
+     &          petab_child(1)-pttab_child(1)+1,
+     &          indmax(2)-indmin(2)+1,
+     &          petab_child(2)-pttab_child(2)+1,
+     &          s_parent(2),s_child(2),ds_parent(2),ds_child(2),2)
+      precomputedone(2) = .TRUE.
+      ENDIF
+
       do k = indmin(nbdim),indmax(nbdim)
 C        
@@ -1490,4 +1597,7 @@
       enddo
       
+      precomputedone(1) = .FALSE.
+      precomputedone(2) = .FALSE.
+      coeffraf = nint ( ds_parent(3) / ds_child(3) )
 
       Call Agrif_Compute_nbdim_interp(s_parent(nbdim),s_child(nbdim),
@@ -1851,5 +1961,6 @@
        ELSEIF (TYPEinterp .EQ. AGRIF_LINEAR) then    
 C
-C         Linear interpolation         
+C         Linear interpolation 
+  
           Call linear1D
      &         (parenttab,childtab,
@@ -1857,4 +1968,12 @@
      &          s_parent,s_child,ds_parent,ds_child)
 C          
+      elseif ( TYPEinterp .EQ. AGRIF_PPM ) then
+
+          Call ppm1D
+     &         (parenttab,childtab,
+     &         indmax-indmin+1,petab_child-pttab_child+1,
+     &         s_parent,s_child,ds_parent,ds_child)
+C
+
         elseif (TYPEinterp .EQ. AGRIF_LAGRANGE) then
 C          
@@ -1906,10 +2025,4 @@
      &         s_parent,s_child,ds_parent,ds_child)
 C              
-      elseif ( TYPEinterp .EQ. AGRIF_PPM ) then
-          Call ppm1D         
-     &         (parenttab,childtab,
-     &         indmax-indmin+1,petab_child-pttab_child+1,
-     &         s_parent,s_child,ds_parent,ds_child)
-C
       endif 
 C
Index: /trunk/AGRIF/AGRIF_FILES/modinterpbasic.F
===================================================================
--- /trunk/AGRIF/AGRIF_FILES/modinterpbasic.F	(revision 778)
+++ /trunk/AGRIF/AGRIF_FILES/modinterpbasic.F	(revision 779)
@@ -37,5 +37,12 @@
 C             
       Real,Dimension(Agrif_MaxRaff) :: tabdiff2, tabdiff3
-      Real,Dimension(:),Allocatable::tabtest4      
+      Real,Dimension(5,Agrif_MaxRaff,3) :: tabppm
+      Real,Dimension(:),Allocatable::tabtest4   
+      Integer, Dimension(:,:), Allocatable :: indparent
+      Integer, Dimension(:,:), Allocatable :: 
+     &           indparentppm,indchildppm
+      Integer, Dimension(:), Allocatable :: 
+     &           indparentppm_1d,indchildppm_1d
+      Real, Dimension(:,:),Allocatable :: coeffparent   
        
       CONTAINS
@@ -140,4 +147,251 @@
 C       
       End Subroutine Linear1d   
+      
+      Subroutine Linear1dprecompute2d(np2, np,nc,
+     &                    s_parent,s_child,ds_parent,ds_child,dir)
+C
+CCC   Description:
+CCC   Subroutine to compute 2D coefficient and index for a linear 1D interpolation
+CCC   on a child grid (vector y) from its parent grid (vector x).
+C
+CC    Method:
+C
+C     Declarations:
+C
+
+C        
+C     Arguments
+      INTEGER             :: np,nc,np2,dir
+      REAL                :: s_parent,s_child,ds_parent,ds_child
+C
+C     Local scalars
+      INTEGER :: i,j,coeffraf,locind_parent_left,inc,inc1,inc2,toto 
+      Integer, Dimension(:,:), Allocatable :: indparent_tmp
+      Real, Dimension(:,:), Allocatable :: coeffparent_tmp
+      REAL    :: ypos,globind_parent_left,globind_parent_right
+      REAL    :: invds, invds2
+      REAL :: ypos2,diff
+C
+C
+
+      coeffraf = nint(ds_parent/ds_child)
+C
+      ypos = s_child
+
+      locind_parent_left = 1 + agrif_int((ypos - s_parent)/ds_parent)
+
+        globind_parent_left = s_parent
+     &                        + (locind_parent_left - 1)*ds_parent
+
+        globind_parent_right = globind_parent_left + ds_parent
+
+C
+      invds = 1./ds_parent
+      invds2 = ds_child/ds_parent
+
+      ypos2 = ypos*invds
+      globind_parent_right=globind_parent_right*invds
+
+      if (.not.allocated(indparent)) then
+      allocate(indparent(nc*np2,3),coeffparent(nc*np2,3))
+      else
+      if (size(indparent,1)<nc*np2) then
+      allocate(coeffparent_tmp(size(indparent,1),size(indparent,2)))
+      allocate(indparent_tmp(size(indparent,1),size(indparent,2)))
+      coeffparent_tmp=coeffparent 
+      indparent_tmp=indparent 
+      deallocate(indparent,coeffparent)
+      allocate(indparent(nc*np2,3),coeffparent(nc*np2,3))
+      coeffparent(1:size(coeffparent_tmp,1),1:size(coeffparent_tmp,2))
+     &  = coeffparent_tmp 
+      indparent(1:size(indparent_tmp,1),1:size(indparent_tmp,2))
+     & = indparent_tmp
+      deallocate(indparent_tmp,coeffparent_tmp)
+      endif
+      endif
+
+      do i = 1,nc-1
+C
+        if (ypos2 > globind_parent_right) then
+           locind_parent_left = locind_parent_left + 1
+           globind_parent_right = globind_parent_right + 1.
+        endif
+        
+        diff=(globind_parent_right - ypos2)
+        indparent(i,dir) = locind_parent_left
+        coeffparent(i,dir) = diff
+        
+        ypos2 = ypos2 + invds2
+C
+      enddo
+C
+      ypos = s_child + (nc-1)*ds_child
+      locind_parent_left = 1 + agrif_int((ypos - s_parent)/ds_parent)
+
+      if (locind_parent_left == np) then
+      indparent(nc,dir) = locind_parent_left-1
+      coeffparent(nc,dir) = 0.
+        else
+C
+          globind_parent_left = s_parent
+     &                        + (locind_parent_left - 1)*ds_parent
+C      
+       indparent(nc,dir) = locind_parent_left
+       
+         coeffparent(nc,dir) = (globind_parent_left + ds_parent - ypos)
+     &       * invds
+      endif                                           
+
+      do i=2, np2
+      inc  =  i*nc
+      inc1 = (i-1)*nc
+      inc2 = (i-2)*nc 
+!CDIR ALTCODE
+      indparent(1+inc1:inc,dir) = indparent(1+inc2:inc1,dir)+np
+!CDIR ALTCODE
+      coeffparent(1+inc1:inc,dir) =coeffparent(1:nc,dir)
+      enddo      
+
+      Return
+C
+C       
+      End Subroutine Linear1dprecompute2d
+
+
+
+      Subroutine Linear1dprecompute(np,nc,
+     &                    s_parent,s_child,ds_parent,ds_child) 
+C
+CCC   Description:
+CCC   Subroutine to compute 1D coefficient and index for a linear 1D interpolation
+CCC   on a child grid (vector y) from its parent grid (vector x).
+C
+CC    Method:
+C
+C     Declarations:
+C
+      
+C        
+C     Arguments
+      INTEGER             :: np,nc 
+      REAL                :: s_parent,s_child,ds_parent,ds_child
+C
+C     Local scalars
+      INTEGER :: i,coeffraf,locind_parent_left
+      REAL    :: ypos,globind_parent_left,globind_parent_right
+      REAL    :: invds, invds2
+      REAL :: ypos2,diff
+C
+C
+
+      coeffraf = nint(ds_parent/ds_child)
+C
+      if (coeffraf == 1) then
+C
+          return
+C
+      endif                          
+C
+      ypos = s_child 
+
+      locind_parent_left = 1 + agrif_int((ypos - s_parent)/ds_parent)
+
+        globind_parent_left = s_parent 
+     &                        + (locind_parent_left - 1)*ds_parent
+
+        globind_parent_right = globind_parent_left + ds_parent
+
+C
+      invds = 1./ds_parent
+      invds2 = ds_child/ds_parent
+      
+      ypos2 = ypos*invds
+      globind_parent_right=globind_parent_right*invds
+      
+      if (.not.allocated(indparent)) then
+      allocate(indparent(nc,1),coeffparent(nc,1))
+      else
+      if (size(indparent)<nc) then
+      deallocate(indparent,coeffparent)
+      allocate(indparent(nc,1),coeffparent(nc,1))
+      endif
+      endif
+      
+      do i = 1,nc-1
+C
+        if (ypos2 > globind_parent_right) then
+           locind_parent_left = locind_parent_left + 1
+           globind_parent_right = globind_parent_right + 1.
+        endif
+        
+        diff=(globind_parent_right - ypos2)
+        indparent(i,1) = locind_parent_left
+        coeffparent(i,1) = diff
+        ypos2 = ypos2 + invds2
+C
+      enddo
+C
+      ypos = s_child + (nc-1)*ds_child
+      locind_parent_left = 1 + agrif_int((ypos - s_parent)/ds_parent)
+      
+      if (locind_parent_left == np) then
+      indparent(nc,1) = locind_parent_left-1
+      coeffparent(nc,1) = 0.
+        else
+C
+          globind_parent_left = s_parent 
+     &                        + (locind_parent_left - 1)*ds_parent  
+C      
+       indparent(nc,1) = locind_parent_left
+       
+         coeffparent(nc,1) = (globind_parent_left + ds_parent - ypos)
+     &       * invds
+      endif                                          
+C           
+      Return
+C
+C       
+      End Subroutine Linear1dprecompute 
+      
+      Subroutine Linear1daftercompute(x,y,np,nc,
+     &                    s_parent,s_child,ds_parent,ds_child,dir) 
+C
+CCC   Description:
+CCC   Subroutine to do a linear 1D interpolation on a child grid (vector y) from
+CCC   its parent grid (vector x) using precomputed coefficient and index.  
+C
+CC    Method:
+C
+C     Declarations:
+C
+      
+C        
+C     Arguments
+      INTEGER             :: np,nc,dir      
+      REAL,INTENT(IN), DIMENSION(np) :: x      
+      REAL,INTENT(OUT), DIMENSION(nc) :: y  
+      REAL                :: s_parent,s_child,ds_parent,ds_child
+C
+C     Local scalars
+      INTEGER :: i,coeffraf,locind_parent_left
+      REAL    :: ypos,globind_parent_left,globind_parent_right
+      REAL    :: invds, invds2
+      REAL :: ypos2,diff
+C
+C
+
+!CDIR ALTCODE
+!CDIR NODEP
+      do i = 1,nc
+C
+       y(i)=coeffparent(i,dir)*x(indparent(i,dir))+
+     &      (1.-coeffparent(i,dir))*x(indparent(i,dir)+1)
+C
+      enddo                                     
+C           
+      Return
+C
+C       
+      End Subroutine Linear1daftercompute            
        
 C
@@ -569,4 +823,6 @@
       If (coeffraf == 1) Then
           locind_parent_left = 1 + nint((s_child - s_parent)/ds_parent)
+!CDIR ALTCODE
+!CDIR SHORTLOOP
           y(1:nc) = x(locind_parent_left:locind_parent_left+nc-1)
           return
@@ -582,5 +838,6 @@
          Allocate(tabtest4(-2*coeffraf:nc+2*coeffraf))
          ENDIF
-      ENDIF  
+      ENDIF
+            
       ypos = s_child  
 C
@@ -595,4 +852,5 @@
 C
       
+!CDIR NOVECTOR
       Do i=1,coeffraf
         tabdiff2(i)=(real(i)-0.5)*invcoeffraf
@@ -602,4 +860,5 @@
       tabdiff3(1) = (1./3.)*a
       a=2.*a
+!CDIR NOVECTOR
       Do i=2,coeffraf
         tabdiff3(i) = tabdiff3(i-1)+(real(i)-1)*a
@@ -621,8 +880,12 @@
 C 
 C
+!CDIR ALTCODE
+!CDIR SHORTLOOP
       Do i = nmin,nmax
          slope(i) = x(i) - x(i-1)
       Enddo
 
+!CDIR ALTCODE
+!CDIR SHORTLOOP
       Do i = nmin+1,nmax-1
          xl(i)= 0.5*(x(i-1)+x(i))
@@ -631,4 +894,6 @@
 C
 C apply parabolic monotonicity
+!CDIR ALTCODE
+!CDIR SHORTLOOP
        Do i = locind_parent_left,locind_parent_last
           delta(i) = xl(i+1) - xl(i)
@@ -644,4 +909,6 @@
         Do iparent = locind_parent_left,locind_parent_last       
              pos=1
+!CDIR ALTCODE
+!CDIR SHORTLOOP
              Do jj = ipos - coeffraf/2+diffmod,ipos + coeffraf/2
 C
@@ -653,11 +920,307 @@
              ipos = ipos + coeffraf
 C
-        End do     
-C
-C
+        End do
+C
+C
+!CDIR ALTCODE
+!CDIR SHORTLOOP
         y(1:nc)=tabtest4(1:nc)
         
       Return
       End Subroutine ppm1d
+      
+
+      Subroutine ppm1dprecompute2d(np2,np,nc,
+     &                    s_parent,s_child,ds_parent,ds_child,dir)
+C
+CCC   Description:
+CCC   Subroutine to compute 2D coefficients and index for a 1D interpolation
+CCC   using piecewise parabolic method  
+CC    Method: 
+C
+C     Declarations:
+C
+      Implicit none
+C        
+C     Arguments
+      Integer             :: np2,np,nc,dir
+C      Real, Dimension(:),Allocatable :: ytemp
+      Real                :: s_parent,s_child,ds_parent,ds_child
+C
+C     Local scalars
+      Integer, Dimension(:,:), Allocatable :: indparent_tmp
+      Integer, Dimension(:,:), Allocatable :: indchild_tmp
+      Integer :: i,coeffraf,locind_parent_left,locind_parent_last
+      Integer :: iparent,ipos,pos,nmin,nmax
+      Real    :: ypos
+      integer :: i1,jj
+      Real :: xpmin,a
+C      
+      Real :: xrmin,xrmax,am3,s2,s1
+      Real, Dimension(np) :: xl,delta,a6,slope
+      INTEGER :: diffmod
+      REAL :: invcoeffraf
+C      
+      coeffraf = nint(ds_parent/ds_child)
+C
+      invcoeffraf = ds_child/ds_parent
+C 
+
+      if (.not.allocated(indparentppm)) then
+      allocate(
+     & indparentppm(np2*nc,3),
+     & indchildppm(np2*nc,3)
+     &        )
+      else
+      if (size(indparentppm,1)<np2*nc) then
+      allocate(
+     &  indparent_tmp(size(indparentppm,1),size(indparentppm,2)),
+     &  indchild_tmp(size(indparentppm,1),size(indparentppm,2)))
+      indparent_tmp = indparentppm 
+      indchild_tmp  = indchildppm 
+      deallocate(indparentppm,indchildppm)
+      allocate(
+     &indparentppm(np2*nc,3),
+     &indchildppm(np2*nc,3)
+     &        )
+      indparentppm(1:size(indparent_tmp,1),1:size(indparent_tmp,2))
+     & = indparent_tmp
+      indchildppm(1:size(indparent_tmp,1),1:size(indparent_tmp,2))
+     & = indchild_tmp
+      deallocate(indparent_tmp,indchild_tmp)
+      endif
+      endif
+
+      if (.not.allocated(indparentppm_1d)) then
+      allocate(indparentppm_1d(-2*coeffraf:nc+2*coeffraf),
+     &         indchildppm_1d(-2*coeffraf:nc+2*coeffraf))
+      else
+      if (size(indparentppm_1d)<nc+4*coeffraf+1) then
+      deallocate(indparentppm_1d,indchildppm_1d)
+      allocate(indparentppm_1d(-2*coeffraf:nc+2*coeffraf),
+     &         indchildppm_1d(-2*coeffraf:nc+2*coeffraf))
+      endif
+      endif
+
+
+      ypos = s_child
+C
+      locind_parent_left = 1 + agrif_int((ypos - s_parent)/ds_parent)
+      locind_parent_last = 1 +
+     &      agrif_ceiling((ypos +(nc - 1)
+     &      *ds_child - s_parent)/ds_parent)
+C
+      xpmin = s_parent + (locind_parent_left-1)*ds_parent
+      i1 = 1+agrif_int((xpmin-s_child)/ds_child)
+C     
+C
+
+      Do i=1,coeffraf
+        tabdiff2(i)=(real(i)-0.5)*invcoeffraf
+      EndDo
+
+      a = invcoeffraf**2
+      tabdiff3(1) = (1./3.)*a
+      a=2.*a
+!CDIR ALTCODE
+      Do i=2,coeffraf
+        tabdiff3(i) = tabdiff3(i-1)+(real(i)-1)*a
+      EndDo
+
+!CDIR ALTCODE
+      Do i=1,coeffraf
+      tabppm(1,i,dir) = 0.08333333333333*
+     &              (-1.+4*tabdiff2(i)-3*tabdiff3(i))
+      tabppm(2,i,dir) = 0.08333333333333*
+     &              (7.-26.*tabdiff2(i)+18.*tabdiff3(i))
+      tabppm(3,i,dir) = 0.08333333333333*
+     &              (7.+30*tabdiff2(i)-30*tabdiff3(i))
+      tabppm(4,i,dir) = 0.08333333333333*
+     &              (-1.-10.*tabdiff2(i)+18.*tabdiff3(i))
+      tabppm(5,i,dir) = 0.08333333333333*
+     &              (2*tabdiff2(i)-3*tabdiff3(i))                   
+      End Do
+C 
+C
+        diffmod = 0
+       IF (mod(coeffraf,2) == 0) diffmod = 1
+C
+        ipos = i1
+C               
+        Do iparent = locind_parent_left,locind_parent_last
+             pos=1
+!CDIR ALTCODE
+             Do jj = ipos - coeffraf/2+diffmod,ipos + coeffraf/2
+         indparentppm_1d(jj) = iparent-2
+         indchildppm_1d(jj) = pos
+               pos = pos+1
+             End do
+             ipos = ipos + coeffraf
+C
+        End do
+     
+       Do i=1,np2
+
+       indparentppm(1+(i-1)*nc:i*nc,dir) = 
+     &             indparentppm_1d(1:nc) + (i-1)*np
+       indchildppm (1+(i-1)*nc:i*nc,dir) = 
+     &             indchildppm_1d (1:nc) + (i-1)*np
+
+       End do
+
+      Return
+      End Subroutine ppm1dprecompute2d
+
+
+      Subroutine ppm1dprecompute(np,nc,
+     &                    s_parent,s_child,ds_parent,ds_child) 
+C
+CCC   Description:
+CCC   Subroutine to compute coefficient and index for  a 1D interpolation 
+CCC   using piecewise parabolic method  
+CC    Method:
+C
+C     Declarations:
+C
+      Implicit none
+C        
+C     Arguments
+      Integer             :: np,nc
+C      Real, Dimension(:),Allocatable :: ytemp
+      Real                :: s_parent,s_child,ds_parent,ds_child
+C
+C     Local scalars
+      Integer :: i,coeffraf,locind_parent_left,locind_parent_last
+      Integer :: iparent,ipos,pos,nmin,nmax
+      Real    :: ypos
+      integer :: i1,jj
+      Real :: xpmin,a
+C      
+      Real :: xrmin,xrmax,am3,s2,s1  
+      Real, Dimension(np) :: xl,delta,a6,slope
+C      Real, Dimension(:),Allocatable  :: diff,diff2,diff3    
+      INTEGER :: diffmod
+      REAL :: invcoeffraf
+C      
+      coeffraf = nint(ds_parent/ds_child)
+C
+      If (coeffraf == 1) Then
+          return
+      End If
+      invcoeffraf = ds_child/ds_parent
+C 
+      
+      if (.not.allocated(indparentppm)) then
+      allocate(indparentppm(-2*coeffraf:nc+2*coeffraf,1),
+     &         indchildppm(-2*coeffraf:nc+2*coeffraf,1))
+      else
+      if (size(indparentppm,1)<nc+4*coeffraf+1) then
+      deallocate(indparentppm,indchildppm)
+      allocate(indparentppm(-2*coeffraf:nc+2*coeffraf,1),
+     &         indchildppm(-2*coeffraf:nc+2*coeffraf,1))
+      endif
+      endif
+            
+      ypos = s_child  
+C
+      locind_parent_left = 1 + agrif_int((ypos - s_parent)/ds_parent) 
+      locind_parent_last = 1 +
+     &      agrif_ceiling((ypos +(nc - 1) 
+     &      *ds_child - s_parent)/ds_parent)  
+C
+      xpmin = s_parent + (locind_parent_left-1)*ds_parent       
+      i1 = 1+agrif_int((xpmin-s_child)/ds_child)        
+C     
+C
+      
+      Do i=1,coeffraf
+        tabdiff2(i)=(real(i)-0.5)*invcoeffraf
+      EndDo
+
+      a = invcoeffraf**2 
+      tabdiff3(1) = (1./3.)*a
+      a=2.*a
+!CDIR ALTCODE
+!!!CDIR SHORTLOOP
+      Do i=2,coeffraf
+        tabdiff3(i) = tabdiff3(i-1)+(real(i)-1)*a
+      EndDo
+      
+!CDIR ALTCODE
+!!!CDIR SHORTLOOP
+      Do i=1,coeffraf
+      tabppm(1,i,1) = 0.08333333333333*(-1.+4*tabdiff2(i)-3*tabdiff3(i))
+      tabppm(2,i,1) = 0.08333333333333*
+     &              (7.-26.*tabdiff2(i)+18.*tabdiff3(i))
+      tabppm(3,i,1) =0.08333333333333*(7.+30*tabdiff2(i)-30*tabdiff3(i))
+      tabppm(4,i,1) = 0.08333333333333*
+     &              (-1.-10.*tabdiff2(i)+18.*tabdiff3(i))     
+      tabppm(5,i,1) = 0.08333333333333*(2*tabdiff2(i)-3*tabdiff3(i)) 
+      End Do   
+C 
+C
+        diffmod = 0
+       IF (mod(coeffraf,2) == 0) diffmod = 1
+C
+        ipos = i1
+C               
+        Do iparent = locind_parent_left,locind_parent_last       
+             pos=1
+!CDIR ALTCODE
+!CDIR SHORTLOOP
+             Do jj = ipos - coeffraf/2+diffmod,ipos + coeffraf/2
+         indparentppm(jj,1) = iparent-2
+         indchildppm(jj,1) = pos
+               pos = pos+1 
+             End do 
+             ipos = ipos + coeffraf
+C
+        End do
+        
+      Return
+      End Subroutine ppm1dprecompute
+      
+      Subroutine ppm1daftercompute(x,y,np,nc,
+     &                    s_parent,s_child,ds_parent,ds_child,dir) 
+C
+CCC   Description:
+CCC   Subroutine to do a 1D interpolation and apply monotonicity constraints
+CCC   using piecewise parabolic method  
+CCC   on a child grid (vector y) from its parent grid (vector x).
+CC    Method: Use precomputed coefficient and index
+C
+C     Declarations:
+C
+      Implicit none
+C        
+C     Arguments
+      Integer             :: np,nc,dir
+      Real, INTENT(IN),Dimension(np) :: x      
+      Real, INTENT(OUT),Dimension(nc) :: y
+C      Real, Dimension(:),Allocatable :: ytemp
+      Real                :: s_parent,s_child,ds_parent,ds_child
+C
+C     Local scalars
+      Integer :: i,coeffraf,locind_parent_left,locind_parent_last
+      Integer :: iparent,ipos,pos,nmin,nmax
+      Real    :: ypos
+      integer :: i1,jj
+      Real :: xpmin,a
+C      
+      Real :: xrmin,xrmax,am3,s2,s1  
+      Real, Dimension(np) :: xl,delta,a6,slope
+      INTEGER :: diffmod
+C      
+C
+      do i=1,nc
+      y(i)=tabppm(1,indchildppm(i,dir),dir)*x(indparentppm(i,dir))+
+     &        tabppm(2,indchildppm(i,dir),dir)*x(indparentppm(i,dir)+1)+
+     &        tabppm(3,indchildppm(i,dir),dir)*x(indparentppm(i,dir)+2)+
+     &        tabppm(4,indchildppm(i,dir),dir)*x(indparentppm(i,dir)+3)+
+     &        tabppm(5,indchildppm(i,dir),dir)*x(indparentppm(i,dir)+4)
+      enddo
+        
+      Return
+      End Subroutine ppm1daftercompute            
       
 C     **************************************************************************  
Index: /trunk/AGRIF/AGRIF_FILES/modmask.F
===================================================================
--- /trunk/AGRIF/AGRIF_FILES/modmask.F	(revision 778)
+++ /trunk/AGRIF/AGRIF_FILES/modmask.F	(revision 779)
@@ -35,5 +35,5 @@
 C
       IMPLICIT NONE
-      Integer, Parameter :: MaxSearch = 5
+      Integer, Parameter :: MaxSearch = 3
 C
       CONTAINS
@@ -103,8 +103,16 @@
          IF (tempP%var%array3(i0,j0,k0)
      &                        == Agrif_SpecialValue) Then
-            Call CalculNewValTempP((/i0,j0,k0/),
-     &                             tempP,parent,
-     &                             ppbtab,ppetab,
-     &                             noraftab,nbdim)
+!------CDIR NEXPAND
+            Call CalculNewValTempP3D((/i0,j0,k0/),
+     &      tempP%var%array3(ppbtab(1),ppbtab(2),ppbtab(3)),
+     &      parent%var%array3(ppbtab(1),ppbtab(2),ppbtab(3)),
+     &                           ppbtab,ppetab,
+     &                           noraftab,MaxSearch,Agrif_SpecialValue)
+     
+c            Call CalculNewValTempP((/i0,j0,k0/),
+c     &                             tempP,parent,
+c     &                             ppbtab,ppetab,
+c     &                             noraftab,nbdim)
+          
          ENDIF
          enddo
@@ -118,8 +126,9 @@
          IF (tempP%var%array4(i0,j0,k0,l0) 
      &                        == Agrif_SpecialValue) Then
-            Call CalculNewValTempP((/i0,j0,k0,l0/),
-     &                             tempP,parent,
-     &                             ppbtab,ppetab,
-     &                             noraftab,nbdim)
+            Call CalculNewValTempP4D((/i0,j0,k0,l0/),
+     &      tempP%var%array4(ppbtab(1),ppbtab(2),ppbtab(3),ppbtab(4)),
+     &      parent%var%array4(ppbtab(1),ppbtab(2),ppbtab(3),ppbtab(4)),
+     &                           ppbtab,ppetab,
+     &                           noraftab,MaxSearch,Agrif_SpecialValue)
          ENDIF
          enddo
@@ -221,5 +230,7 @@
       Valmax = min(Valmax,MaxSearch)
 C
+!CDIR NOVECTOR
       imin = indic
+!CDIR NOVECTOR
       imax = indic
 C
@@ -238,4 +249,5 @@
                   if (indic(iii).GT.ppbtab(iii)) then
 
+!CDIR NOVECTOR
                      idecal = indic
                      idecal(iii) = idecal(iii)-1
@@ -286,4 +298,6 @@
             SELECT CASE(nbdim)
             CASE (1)
+!CDIR ALTCODE
+!CDIR SHORTLOOP
                do ii = imin(1),imax(1)
                     ValParent = parent%var%array1(ii)
@@ -296,4 +310,6 @@
             CASE (2)
                do jj = imin(2),imax(2)
+!CDIR ALTCODE
+!CDIR SHORTLOOP
                do ii = imin(1),imax(1)
                     ValParent = parent%var%array2(ii,jj)
@@ -308,4 +324,6 @@
                do kk = imin(3),imax(3)
                do jj = imin(2),imax(2)
+!CDIR ALTCODE
+!CDIR SHORTLOOP
                do ii = imin(1),imax(1)
                     ValParent = parent%var%array3(ii,jj,kk)
@@ -322,4 +340,6 @@
                do kk = imin(3),imax(3)
                do jj = imin(2),imax(2)
+!CDIR ALTCODE
+!CDIR SHORTLOOP
                do ii = imin(1),imax(1)
                     ValParent = parent%var%array4(ii,jj,kk,ll)
@@ -338,4 +358,6 @@
                do kk = imin(3),imax(3)
                do jj = imin(2),imax(2)
+!CDIR ALTCODE
+!CDIR SHORTLOOP
                do ii = imin(1),imax(1)
                     ValParent = parent%var%array5(ii,jj,kk,ll,mm)
@@ -356,4 +378,6 @@
                do kk = imin(3),imax(3)
                do jj = imin(2),imax(2)
+!CDIR ALTCODE
+!CDIR SHORTLOOP
                do ii = imin(1),imax(1)
                     ValParent = parent%var%array6(ii,jj,kk,ll,mm,nn)
@@ -412,3 +436,254 @@
 C
 C
-      End Module Agrif_Mask       
+      End Module Agrif_Mask     
+      
+      Subroutine CalculNewValTempP3D(indic,
+     &               tempP,parent,ppbtab,
+     &               ppetab,noraftab,MaxSearch,Agrif_SpecialValue)
+C
+CCC   Description:
+CCC   Subroutine called in the procedure Agrif_InterpnD to recalculate the value 
+CCC   of the parent grid variable when this one is equal to Agrif_SpecialValue. 
+C
+C     Declarations:
+C
+       
+C
+C     Arrays arguments      
+      INTEGER, PARAMETER :: nbdim = 3
+      INTEGER,DIMENSION(nbdim) :: indic
+      LOGICAL,DIMENSION(nbdim) :: noraftab
+      INTEGER,DIMENSION(nbdim) :: ppbtab,ppetab
+      REAL :: Agrif_SpecialValue
+C
+C     Pointer argument
+      REAL,DIMENSION(ppbtab(1):ppetab(1),ppbtab(2):ppetab(2),
+     &    ppbtab(3):ppetab(3)) :: tempP, parent  ! Part of the parent grid used for 
+                                      ! the interpolation of the child grid                     
+C
+C     Local scalar
+      INTEGER                  :: i,ii,iii,jj,kk,ll,mm,nn 
+      INTEGER,DIMENSION(nbdim) :: imin,imax,idecal
+      INTEGER                  :: Nbvals
+      REAL                     :: Res
+      REAL                     :: ValParent
+      INTEGER                  :: ValMax
+      INTEGER :: MaxSearch
+      LOGICAL :: Existunmasked
+C      
+C     Local arrays      
+C
+      ValMax = 1
+!CDIR NOVECTOR
+      do iii = 1 , nbdim
+         IF (.NOT.noraftab(iii)) THEN
+         ValMax = max(ValMax,ppetab(iii)-indic(iii))
+         ValMax = max(ValMax,indic(iii)-ppbtab(iii))
+         ENDIF
+      enddo
+C
+      Valmax = min(Valmax,MaxSearch)
+C
+!CDIR NOVECTOR
+      imin = indic
+!CDIR NOVECTOR
+      imax = indic
+      
+!CDIR NOVECTOR
+         idecal = indic	 
+C
+         i = Valmax
+
+            do iii = 1 , nbdim
+               if (.NOT.noraftab(iii)) then
+                  imin(iii) = max(indic(iii) - i,ppbtab(iii))
+                  imax(iii) = min(indic(iii) + i,ppetab(iii))
+         
+                  if (indic(iii).GT.ppbtab(iii)) then
+
+                     idecal(iii) = idecal(iii)-1
+
+                        if (tempP(idecal(1),idecal(2),idecal(3)) 
+     &                               == Agrif_SpecialValue) then
+                           imin(iii) = imax(iii)
+                        endif
+
+                     idecal(iii) = idecal(iii)+1
+                  endif
+
+               endif            
+            enddo
+C
+            Existunmasked = .FALSE.
+C
+               do kk = imin(3),imax(3)
+               do jj = imin(2),imax(2)
+!CDIR NOVECTOR
+               do ii = imin(1),imax(1)
+                    if ( parent(ii,jj,kk) .NE. Agrif_SpecialValue) then
+                        Existunmasked = .TRUE.
+                       exit
+                    endif 
+                        enddo
+                  enddo  
+               enddo
+C
+C
+          if (.Not.Existunmasked) return
+C
+         i = 1
+C
+         do While(i <= ValMax)
+C
+
+            do iii = 1 , nbdim
+               if (.NOT.noraftab(iii)) then
+                  imin(iii) = max(indic(iii) - i,ppbtab(iii))
+                  imax(iii) = min(indic(iii) + i,ppetab(iii))
+               endif            
+            enddo
+C
+            Res = 0.
+            Nbvals = 0
+C
+               do kk = imin(3),imax(3)
+               do jj = imin(2),imax(2)
+!CDIR NOVECTOR
+               do ii = imin(1),imax(1)
+                    ValParent = parent(ii,jj,kk)
+                    if ( ValParent .NE. Agrif_SpecialValue) then
+                        Res = Res + ValParent
+                        Nbvals = Nbvals + 1
+                    endif 
+                        enddo
+                  enddo  
+               enddo
+C
+C
+            
+            if (Nbvals.GT.0) then
+              tempP(indic(1),indic(2),indic(3)) = Res/Nbvals
+              exit
+            else
+               i = i + 1                      
+            endif
+          enddo            
+C     
+      End Subroutine CalculNewValTempP3D        
+
+      Subroutine CalculNewValTempP4D(indic,
+     &               tempP,parent,ppbtab,
+     &               ppetab,noraftab,MaxSearch,Agrif_SpecialValue)
+C
+CCC   Description:
+CCC   Subroutine called in the procedure Agrif_InterpnD to recalculate the value 
+CCC   of the parent grid variable when this one is equal to Agrif_SpecialValue. 
+C
+C     Declarations:
+C
+       
+C
+C     Arrays arguments      
+      INTEGER, PARAMETER :: nbdim = 4
+      INTEGER,DIMENSION(nbdim) :: indic
+      LOGICAL,DIMENSION(nbdim) :: noraftab
+      INTEGER,DIMENSION(nbdim) :: ppbtab,ppetab
+      INTEGER :: MaxSearch      
+      REAL :: Agrif_SpecialValue      
+C
+C     Pointer argument
+      REAL,DIMENSION(ppbtab(1):ppetab(1),ppbtab(2):ppetab(2),
+     &    ppbtab(3):ppetab(3),
+     &    ppbtab(4):ppetab(4)) :: tempP, parent  ! Part of the parent grid used for 
+                                      ! the interpolation of the child grid                     
+C
+C     Local scalar
+      INTEGER                  :: i,ii,iii,jj,kk,ll,mm,nn 
+      INTEGER,DIMENSION(nbdim) :: imin,imax,idecal
+      INTEGER                  :: Nbvals
+      REAL                     :: Res
+      REAL                     :: ValParent
+      INTEGER                  :: ValMax
+      LOGICAL                  :: firsttest
+C      
+C     Local arrays      
+C
+      ValMax = 1
+      do iii = 1 , nbdim
+         IF (.NOT.noraftab(iii)) THEN
+         ValMax = max(ValMax,ppetab(iii)-indic(iii))
+         ValMax = max(ValMax,indic(iii)-ppbtab(iii))
+         ENDIF
+      enddo
+C
+      Valmax = min(Valmax,MaxSearch)
+C
+      imin = indic
+      imax = indic
+C
+         i = 1
+         firsttest = .TRUE.
+         idecal = indic	 
+	 
+C
+         do While(i <= ValMax)
+C
+         IF ((i == 1).AND.(firsttest)) i = Valmax
+
+            do iii = 1 , nbdim
+               if (.NOT.noraftab(iii)) then
+                  imin(iii) = max(indic(iii) - i,ppbtab(iii))
+                  imax(iii) = min(indic(iii) + i,ppetab(iii))
+                  if (firsttest) then               
+                  if (indic(iii).GT.ppbtab(iii)) then
+
+
+                     idecal(iii) = idecal(iii)-1
+		     
+                     if (tempP(idecal(1),idecal(2),idecal(3),idecal(4))
+     &                               == Agrif_SpecialValue) then
+                           imin(iii) = imax(iii)
+                     endif
+		     
+		     idecal(iii) = idecal(iii)+1
+                  endif
+                  endif
+               endif            
+            enddo
+C
+            Res = 0.
+            Nbvals = 0
+C
+               do ll = imin(4),imax(4)
+               do kk = imin(3),imax(3)
+               do jj = imin(2),imax(2)
+               do ii = imin(1),imax(1)
+                    ValParent = parent(ii,jj,kk,ll)
+                    if ( ValParent .NE. Agrif_SpecialValue) then
+                        Res = Res + ValParent
+                        Nbvals = Nbvals + 1
+                    endif 
+                              enddo
+                        enddo
+                  enddo  
+               enddo
+C
+C
+            
+            if (Nbvals.GT.0) then
+              if (firsttest) then
+                   firsttest = .FALSE.
+                   i=1
+                   cycle
+              endif
+	      
+              tempP(indic(1),indic(2),indic(3),indic(4)) = Res/Nbvals
+	    
+              exit
+            else
+               if (firsttest) exit
+               i = i + 1                      
+            endif
+          enddo            
+C     
+      End Subroutine CalculNewValTempP4D        
Index: /trunk/AGRIF/AGRIF_FILES/modmpp.F
===================================================================
--- /trunk/AGRIF/AGRIF_FILES/modmpp.F	(revision 778)
+++ /trunk/AGRIF/AGRIF_FILES/modmpp.F	(revision 779)
@@ -146,4 +146,5 @@
         sendtoproc(k) = .true.
 C
+!CDIR SHORTLOOP
         do i = 1,nbdim
 C
@@ -228,4 +229,5 @@
                 datasize = 1
 C
+!CDIR SHORTLOOP
                 do i = 1,nbdim
 C
@@ -305,4 +307,5 @@
                 datasize = 1
 C
+!CDIR SHORTLOOP
                 do i = 1,nbdim
 C
Index: /trunk/AGRIF/AGRIF_FILES/modupdate.F
===================================================================
--- /trunk/AGRIF/AGRIF_FILES/modupdate.F	(revision 778)
+++ /trunk/AGRIF/AGRIF_FILES/modupdate.F	(revision 779)
@@ -42,4 +42,5 @@
 C
       IMPLICIT NONE
+      logical, private :: precomputedone(7) = .FALSE.      
 C      
       CONTAINS
@@ -1301,4 +1302,5 @@
 
       if ( nbdim .EQ. 1 ) then
+         tempP%var%array1 = 0.
          Call Agrif_Update_1D_recursive(TypeUpdate,
      &           tempP%var%array1,tempCextend%var%array1,
@@ -1822,9 +1824,47 @@
       REAL, DIMENSION(indmin(1):indmax(1),
      &                 pttab_child(2):petab_child(2)) :: tabtemp
+      REAL, DIMENSION(indmin(2):indmax(2),
+     &                indmin(1):indmax(1))           :: tempP_trsp
+      REAL, DIMENSION(pttab_child(2):petab_child(2),
+     &                 indmin(1):indmax(1)) :: tabtemp_trsp
       INTEGER :: i,j
       INTEGER :: coeffraf,locind_child_left
 C
+      tabtemp = 0.
+
+
+      coeffraf = nint ( ds_parent(1) / ds_child(1) ) 
+      IF((TypeUpdate(1) == AGRIF_Update_average) 
+     &                 .AND. (coeffraf /= 1 ))THEN
+!---CDIR NEXPAND
+         IF(.NOT. precomputedone(1) ) Call average1Dprecompute2D
+     &       (petab_child(2)-pttab_child(2)+1,
+     &        indmax(1)-indmin(1)+1,
+     &        petab_child(1)-pttab_child(1)+1,
+     &        s_parent(1),s_child(1),ds_parent(1),ds_child(1),1)
+!---CDIR NEXPAND
+          Call average1Daftercompute
+     &       (  tabtemp, tempC,
+     &          size(tabtemp), size(tempC),
+     &          s_parent(1),s_child(1),ds_parent(1),ds_child(1),1)
+ 
+      ELSEIF((TypeUpdate(1) == AGRIF_Update_copy)
+     &                 .AND. (coeffraf /= 1 ))THEN
+!---CDIR NEXPAND
+         IF(.NOT. precomputedone(1) ) Call copy1Dprecompute2D
+     &       (petab_child(2)-pttab_child(2)+1,
+     &        indmax(1)-indmin(1)+1,
+     &        petab_child(1)-pttab_child(1)+1,
+     &        s_parent(1),s_child(1),ds_parent(1),ds_child(1),1)
+!---CDIR NEXPAND
+          Call copy1Daftercompute
+     &       (  tabtemp, tempC,
+     &          size(tabtemp), size(tempC),
+     &          s_parent(1),s_child(1),ds_parent(1),ds_child(1),1)
+
+      ELSE
       do j = pttab_child(nbdim),petab_child(nbdim)
 C
+!---CDIR NEXPAND
         Call Agrif_Update_1D_recursive(TypeUpdate(1:nbdim-1),    
      &         tabtemp(:,j),
@@ -1836,13 +1876,51 @@
 C
       enddo
-      
+      ENDIF
+      tabtemp_trsp = TRANSPOSE(tabtemp) 
+      coeffraf = nint(ds_parent(nbdim)/ds_child(nbdim))
+    
+!---CDIR NEXPAND
       Call Agrif_Compute_nbdim_update(s_parent(nbdim),s_child(nbdim),
      &  ds_parent(nbdim),ds_child(nbdim),coeffraf,locind_child_left)
 C
+      
+      tempP_trsp = 0. 
+
+      IF((TypeUpdate(2) == AGRIF_Update_average)
+     &                .AND. (coeffraf /= 1 ))THEN
+!---CDIR NEXPAND
+         IF(.NOT. precomputedone(2) ) Call average1Dprecompute2D
+     &      ( indmax(1)-indmin(1)+1,
+     &        indmax(2)-indmin(2)+1,
+     &        petab_child(2)-pttab_child(2)+1,
+     &        s_parent(2),s_child(2),ds_parent(2),ds_child(2),2)
+!---CDIR NEXPAND
+          Call average1Daftercompute
+     &       (  tempP_trsp, tabtemp_trsp,
+     &          size(tempP_trsp), size(tabtemp_trsp),
+     &          s_parent(2),s_child(2),ds_parent(2),ds_child(2),2)
+
+      ELSEIF((TypeUpdate(2) == AGRIF_Update_copy)
+     &                .AND. (coeffraf /= 1 ))THEN
+!---CDIR NEXPAND
+         IF(.NOT. precomputedone(2) ) Call copy1Dprecompute2D
+     &      ( indmax(1)-indmin(1)+1,
+     &        indmax(2)-indmin(2)+1,
+     &        petab_child(2)-pttab_child(2)+1,
+     &        s_parent(2),s_child(2),ds_parent(2),ds_child(2),2)
+!---CDIR NEXPAND
+          Call copy1Daftercompute
+     &       (  tempP_trsp, tabtemp_trsp,
+     &          size(tempP_trsp), size(tabtemp_trsp),
+     &          s_parent(2),s_child(2),ds_parent(2),ds_child(2),2)
+
+      ELSE
+
       do i = indmin(1),indmax(1)
 C
+!---CDIR NEXPAND
         Call Agrif_UpdateBase(TypeUpdate(2),
-     &           tempP(i,:),
-     &          tabtemp(i,:),
+     &           tempP_trsp(indmin(nbdim):indmax(nbdim),i),
+     &          tabtemp_trsp(pttab_child(nbdim):petab_child(nbdim),i),
      &           indmin(nbdim),indmax(nbdim),
      &           pttab_child(nbdim),petab_child(nbdim),
@@ -1852,4 +1930,8 @@
 C        
       enddo
+      
+      ENDIF
+
+      tempP = TRANSPOSE(tempP_trsp)
 C
       Return
@@ -1902,4 +1984,47 @@
 C
 C
+      coeffraf = nint ( ds_parent(1) / ds_child(1) )
+      IF((TypeUpdate(1) == AGRIF_Update_average)
+     &                 .AND. (coeffraf /= 1 ))THEN
+!---CDIR NEXPAND
+         Call average1Dprecompute2D
+     &       (petab_child(2)-pttab_child(2)+1,
+     &        indmax(1)-indmin(1)+1,
+     &        petab_child(1)-pttab_child(1)+1,
+     &        s_parent(1),s_child(1),ds_parent(1),ds_child(1),1)
+      precomputedone(1) = .TRUE.
+      ELSEIF((TypeUpdate(1) == AGRIF_Update_copy)
+     &                 .AND. (coeffraf /= 1 ))THEN
+!---CDIR NEXPAND
+         Call copy1Dprecompute2D
+     &       (petab_child(2)-pttab_child(2)+1,
+     &        indmax(1)-indmin(1)+1,
+     &        petab_child(1)-pttab_child(1)+1,
+     &        s_parent(1),s_child(1),ds_parent(1),ds_child(1),1)
+      precomputedone(1) = .TRUE.
+      ENDIF
+
+      coeffraf = nint ( ds_parent(2) / ds_child(2) )
+      IF((TypeUpdate(2) == AGRIF_Update_average)
+     &                .AND. (coeffraf /= 1 ))THEN
+!---CDIR NEXPAND
+         Call average1Dprecompute2D
+     &      ( indmax(1)-indmin(1)+1,
+     &        indmax(2)-indmin(2)+1,
+     &        petab_child(2)-pttab_child(2)+1,
+     &        s_parent(2),s_child(2),ds_parent(2),ds_child(2),2)
+      precomputedone(2) = .TRUE.
+      ELSEIF((TypeUpdate(2) == AGRIF_Update_copy)
+     &                .AND. (coeffraf /= 1 ))THEN
+!---CDIR NEXPAND
+         Call copy1Dprecompute2D
+     &      ( indmax(1)-indmin(1)+1,
+     &        indmax(2)-indmin(2)+1,
+     &        petab_child(2)-pttab_child(2)+1,
+     &        s_parent(2),s_child(2),ds_parent(2),ds_child(2),2)
+      precomputedone(2) = .TRUE.
+      ENDIF
+
+
       do k = pttab_child(nbdim),petab_child(nbdim)
 C
@@ -1914,4 +2039,8 @@
       enddo
 C
+      precomputedone(1) = .FALSE.
+      precomputedone(2) = .FALSE.
+      coeffraf = nint ( ds_parent(3) / ds_child(3) )
+
       Call Agrif_Compute_nbdim_update(s_parent(nbdim),s_child(nbdim),
      &  ds_parent(nbdim),ds_child(nbdim),coeffraf,locind_child_left)
@@ -1930,4 +2059,5 @@
               
       ELSE
+      tempP = 0.
 C
       do j = indmin(2),indmax(2)
@@ -2021,4 +2151,6 @@
      &  ds_parent(nbdim),ds_child(nbdim),coeffraf,locind_child_left)
 C
+      tempP = 0.
+
       do k = indmin(3),indmax(3)
 C
@@ -2120,4 +2252,5 @@
       Call Agrif_Compute_nbdim_update(s_parent(nbdim),s_child(nbdim),
      &  ds_parent(nbdim),ds_child(nbdim),coeffraf,locind_child_left)
+      tempP = 0.
 C
       do l = indmin(4),indmax(4)
@@ -2232,4 +2365,6 @@
      &  ds_parent(nbdim),ds_child(nbdim),coeffraf,locind_child_left)
 C
+      tempP = 0.
+
       do m = indmin(5),indmax(5)
       do l = indmin(4),indmax(4)
@@ -2307,9 +2442,9 @@
 C
         elseif (TypeUpdate == AGRIF_Update_average) then
-C             
+C           
           Call average1D
      &       (parenttab,childtab,
      &          indmax-indmin+1,petab_child-pttab_child+1,
-     &          s_parent,s_child,ds_parent,ds_child)    
+     &          s_parent,s_child,ds_parent,ds_child) 
 C
         elseif (TypeUpdate == AGRIF_Update_full_weighting) then
@@ -2421,5 +2556,4 @@
 #endif
            
-      End Module Agrif_Update
 
 C     **************************************************************************
@@ -2443,5 +2577,4 @@
 C
 C     Arguments
-      USE Agrif_Update
       INTEGER                   :: nbdim
       INTEGER, DIMENSION(nbdim) :: TypeUpdate            ! TYPE of update (copy or average)
@@ -2466,4 +2599,5 @@
      &                  ds_parent(nbdim),ds_child(nbdim),
      &                  coeffraf,locind_child_left)
+
 C
       Return
@@ -2471,2 +2605,4 @@
 C
       End Subroutine Agrif_Update_1D_recursive      
+
+      End Module Agrif_Update
Index: /trunk/AGRIF/AGRIF_FILES/modupdatebasic.F
===================================================================
--- /trunk/AGRIF/AGRIF_FILES/modupdatebasic.F	(revision 778)
+++ /trunk/AGRIF/AGRIF_FILES/modupdatebasic.F	(revision 779)
@@ -37,4 +37,7 @@
       
       IMPLICIT NONE
+      
+      integer,dimension(:,:),allocatable :: indchildcopy
+      integer,dimension(:,:),allocatable :: indchildaverage
 C             
 
@@ -78,4 +81,5 @@
           locind_child_left = 1 + nint((s_parent - s_child)/ds_child)
 C        
+!CDIR ALTCODE
           x(1:np) = y(locind_child_left:locind_child_left+np-1)
 C
@@ -87,4 +91,5 @@
       locind_child_left = 1 + nint((s_parent - s_child)/ds_child)
       
+!CDIR ALTCODE
       do i = 1,np
 C      
@@ -100,6 +105,172 @@
 C
       End Subroutine agrif_copy1d
-C
-C
+
+C     **************************************************************************  
+CCC   Subroutine Copy1dprecompute  
+C     ************************************************************************** 
+C
+      Subroutine copy1dprecompute2d(nc2,np,nc,
+     &                  s_parent,s_child,ds_parent,ds_child,dir)
+C
+CCC   Description:
+CCC   Subroutine to precompute index for a copy on a parent grid (vector x) from
+CCC   its child grid (vector y).
+C
+CC    Method:
+C
+C     Declarations:
+C
+
+C        
+C     Arguments
+      INTEGER             :: nc2,np,nc
+      INTEGER             :: dir   
+      REAL                :: s_parent,s_child
+      REAL                :: ds_parent,ds_child
+C
+C     Local variables
+      INTEGER,DIMENSION(:,:),ALLOCATABLE :: indchildcopy_tmp
+      INTEGER :: i,locind_child_left,coeffraf
+C 
+C
+      coeffraf = nint(ds_parent/ds_child)
+C
+      locind_child_left = 1 + nint((s_parent - s_child)/ds_child)
+
+      if (.not.allocated(indchildcopy)) then
+      allocate(indchildcopy(np*nc2,3))
+      else
+      if (size(indchildcopy,1)<np*nc2) then
+      allocate( indchildcopy_tmp(
+     &         size(indchildcopy,1),size(indchildcopy,2)))
+      indchildcopy_tmp = indchildcopy
+      deallocate(indchildcopy)
+      allocate(indchildcopy(np*nc2,3))
+      indchildcopy(1:size(indchildcopy_tmp,1),
+     &              1:size(indchildcopy_tmp,2)) = indchildcopy_tmp
+      deallocate(indchildcopy_tmp)
+      endif
+      endif
+
+
+      do i = 1,np
+C      
+         indchildcopy(i,dir) = locind_child_left
+         locind_child_left = locind_child_left + coeffraf
+C         
+      enddo
+
+      do i =2, nc2
+        indchildcopy(1+(i-1)*np:i*np,dir)=
+     &       indchildcopy(1+(i-2)*np:(i-1)*np,dir) + nc 
+      enddo
+
+C
+      Return
+C
+C
+      End Subroutine copy1dprecompute2d
+
+C
+C
+C     **************************************************************************  
+CCC   Subroutine Copy1dprecompute  
+C     ************************************************************************** 
+C
+      Subroutine copy1dprecompute(np,nc,
+     &                  s_parent,s_child,ds_parent,ds_child,dir)
+C
+CCC   Description:
+CCC   Subroutine to precompute index for a copy on a parent grid (vector x) from
+CCC   its child grid (vector y).
+C
+CC    Method:
+C
+C     Declarations:
+C
+
+C        
+C     Arguments
+      INTEGER             :: np,nc
+      INTEGER             :: dir
+      REAL                :: s_parent,s_child
+      REAL                :: ds_parent,ds_child
+C
+C     Local variables
+      INTEGER :: i,locind_child_left,coeffraf
+C 
+C
+      coeffraf = nint(ds_parent/ds_child)
+C
+      locind_child_left = 1 + nint((s_parent - s_child)/ds_child)
+
+      if (.not.allocated(indchildcopy)) then
+      allocate(indchildcopy(np,1))
+      else
+      if (size(indchildcopy)<np) then
+      deallocate(indchildcopy)
+      allocate(indchildcopy(np,1))
+      endif
+      endif
+
+!CDIR ALTCODE
+      do i = 1,np
+C      
+         indchildcopy(i,1) = locind_child_left
+         locind_child_left = locind_child_left + coeffraf
+C         
+      enddo
+
+C
+      Return
+C
+C
+      End Subroutine copy1dprecompute
+C
+C
+C     **************************************************************************  
+CCC   Subroutine Copy1daftercompute  
+C     ************************************************************************** 
+C
+      Subroutine copy1daftercompute(x,y,np,nc,
+     &                  s_parent,s_child,ds_parent,ds_child,dir)
+C
+CCC   Description:
+CCC   Subroutine to do a copy on a parent grid (vector x) from its child grid 
+CCC   (vector y) using precomputed index.  
+C
+CC    Method:
+C
+C     Declarations:
+C
+  
+C        
+C     Arguments
+      INTEGER             :: np,nc
+      INTEGER             :: dir
+      REAL, DIMENSION(np) :: x
+      REAL, DIMENSION(nc) :: y
+      REAL                :: s_parent,s_child
+      REAL                :: ds_parent,ds_child
+C
+C     Local variables
+      INTEGER :: i
+C 
+C
+      
+!CDIR ALTCODE
+      do i = 1,np
+C      
+         x(i) = y(indchildcopy(i,dir))
+C
+      enddo
+  
+C
+      Return
+C
+C
+      End Subroutine copy1daftercompute
+
+
 C
 C     **************************************************************************  
@@ -123,5 +294,5 @@
 C     Local variables
       INTEGER :: i,locind_child_left,coeffraf,ii
-      REAL    :: xpos 
+      REAL    :: xpos, invcoeffraf 
       INTEGER :: nbnonnuls
       INTEGER :: diffmod
@@ -129,4 +300,5 @@
 C
       coeffraf = nint(ds_parent/ds_child)
+      invcoeffraf = 1./coeffraf
 C
       if (coeffraf == 1) then
@@ -149,49 +321,224 @@
 
         locind_child_left = 1 + agrif_int((xpos - s_child)/ds_child)
-      
-      do i = 1,np
-C
-
-C
-c        if ((locind_child_left-1 < 1) 
-c     &      .OR. (locind_child_left+1 > nc)) then 
-C
-c            x(i) = y(locind_child_left)                
-C
-c          else 
+C
+      IF (Agrif_UseSpecialValueInUpdate) THEN
+
+      do i = 1,np
           nbnonnuls = 0
+!CDIR NOVECTOR
           Do ii = -coeffraf/2+locind_child_left+diffmod,
      &                coeffraf/2+locind_child_left
 C 
-            IF (Agrif_UseSpecialValueInUpdate) THEN
-            IF (y(ii) .NE. Agrif_SpecialValueFineGrid) THEN
+               IF (y(ii) .NE. Agrif_SpecialValueFineGrid) THEN
                nbnonnuls = nbnonnuls + 1
                x(i) = x(i) + y(ii)
-            ENDIF
-            ELSE
+               ENDIF
+          End Do
+               IF (nbnonnuls .NE. 0) THEN
+                  x(i) = x(i)/nbnonnuls
+               ELSE
+                  x(i) = Agrif_SpecialValueFineGrid
+               ENDIF
+        locind_child_left = locind_child_left + coeffraf
+      enddo
+
+      ELSE 
+
+!CDIR ALTCODE
+      do i = 1,np
+!CDIR NOVECTOR
+          Do ii = -coeffraf/2+locind_child_left+diffmod,
+     &                coeffraf/2+locind_child_left
+C 
                x(i) = x(i) + y(ii)
-            ENDIF
           End Do
-            IF (Agrif_UseSpecialValueInUpdate) THEN
-                 IF (nbnonnuls .NE. 0) THEN
-                    x(i) = x(i)/nbnonnuls
-                 ELSE
-                    x(i) = Agrif_SpecialValueFineGrid
-                 ENDIF
-            ELSE
-                 x(i) = x(i)/coeffraf
-            ENDIF
-C
-c       endif
-C
-c        xpos = xpos + ds_parent
+                 x(i) = x(i)*invcoeffraf
         locind_child_left = locind_child_left + coeffraf
-C
-      enddo
-C
+      enddo
+
+      ENDIF
+
       Return 
 C            
 C     
       End Subroutine average1d
+      
+      Subroutine average1dprecompute2d(nc2,np,nc,
+     &                     s_parent,s_child,ds_parent,ds_child,dir)
+C
+CCC   Description:
+CCC   Subroutine to do an update by average on a parent grid (vector x)from its 
+CCC   child grid (vector y).
+C
+C     Arguments
+      INTEGER             :: nc2,np,nc,dir
+      REAL                :: s_parent,s_child
+      REAL                :: ds_parent,ds_child
+C
+C     Local variables
+      INTEGER :: i,locind_child_left,coeffraf,ii
+      INTEGER,DIMENSION(:,:),ALLOCATABLE :: indchildaverage_tmp
+      REAL    :: xpos, invcoeffraf
+      INTEGER :: nbnonnuls
+      INTEGER :: diffmod
+C 
+C
+      coeffraf = nint(ds_parent/ds_child)
+C
+      xpos = s_parent
+C
+      diffmod = 0
+
+      IF ( mod(coeffraf,2) == 0 ) diffmod = 1
+
+      locind_child_left = 1 + agrif_int((xpos - s_child)/ds_child)
+
+      if (.not.allocated(indchildaverage)) then
+      allocate(indchildaverage(np*nc2,3))
+      else
+      if (size(indchildaverage,1)<np*nc2) then
+      allocate( indchildaverage_tmp(
+     &         size(indchildaverage,1),size(indchildaverage,2)))
+      indchildaverage_tmp = indchildaverage
+      deallocate(indchildaverage)
+      allocate(indchildaverage(np*nc2,3))
+      indchildaverage(1:size(indchildaverage_tmp,1),
+     &              1:size(indchildaverage_tmp,2)) = indchildaverage_tmp
+      deallocate(indchildaverage_tmp)
+      endif
+      endif
+
+      do i = 1,np
+        indchildaverage(i,dir)= -coeffraf/2+locind_child_left
+     &                                    +diffmod
+        locind_child_left = locind_child_left + coeffraf
+      enddo
+      do i = 2,nc2
+        indchildaverage(1+(i-1)*np:i*np,dir)= 
+     &                 indchildaverage(1+(i-2)*np:(i-1)*np,dir) + nc 
+      enddo
+
+      Return 
+C            
+C     
+      End Subroutine average1dprecompute2d
+      
+
+      Subroutine average1dprecompute(np,nc,
+     &                     s_parent,s_child,ds_parent,ds_child) 
+C
+CCC   Description:
+CCC   Subroutine to do an update by average on a parent grid (vector x)from its 
+CCC   child grid (vector y).
+C
+C     Arguments
+      INTEGER             :: np,nc 
+      REAL                :: s_parent,s_child
+      REAL                :: ds_parent,ds_child
+C
+C     Local variables
+      INTEGER :: i,locind_child_left,coeffraf,ii
+      REAL    :: xpos, invcoeffraf 
+      INTEGER :: nbnonnuls
+      INTEGER :: diffmod
+C 
+C
+      coeffraf = nint(ds_parent/ds_child)
+
+C
+      if (coeffraf == 1) then
+C
+          return
+C
+      endif
+            
+C
+      xpos = s_parent
+C
+      diffmod = 0
+      
+      IF ( mod(coeffraf,2) == 0 ) diffmod = 1
+
+        locind_child_left = 1 + agrif_int((xpos - s_child)/ds_child)
+
+      if (.not.allocated(indchildaverage)) then
+      allocate(indchildaverage(np,1))
+      else
+      if (size(indchildaverage,1)<np) then
+      deallocate(indchildaverage)
+      allocate(indchildaverage(np,1))
+      endif
+      endif
+      
+!CDIR ALTCODE
+      do i = 1,np
+C 
+        indchildaverage(i,1)=-coeffraf/2+locind_child_left+diffmod
+        locind_child_left = locind_child_left + coeffraf
+      enddo
+
+      Return 
+C            
+C     
+      End Subroutine average1dprecompute
+      
+      Subroutine average1daftercompute(x,y,np,nc,
+     &                     s_parent,s_child,ds_parent,ds_child,dir) 
+C
+CCC   Description:
+CCC   Subroutine to do an update by average on a parent grid (vector x)from its 
+CCC   child grid (vector y).
+C
+C     Arguments
+      INTEGER             :: np,nc,dir      
+      REAL, DIMENSION(np) :: x      
+      REAL, DIMENSION(nc) :: y  
+      REAL                :: s_parent,s_child
+      REAL                :: ds_parent,ds_child
+C
+C     Local variables
+      INTEGER :: i,locind_child_left,coeffraf,ii,j
+      REAL    :: xpos, invcoeffraf 
+      REAL, PARAMETER :: one_third=1./3.
+      INTEGER, DIMENSION(np) :: nbnonnuls
+      INTEGER :: diffmod
+      REAL, DIMENSION(0:5) :: invcoeff=(/1.,1.,0.5,one_third,0.25,0.2/)
+C 
+C
+      coeffraf = nint(ds_parent/ds_child)
+      invcoeffraf = 1./coeffraf      
+C      
+C
+      IF (Agrif_UseSpecialValueInUpdate) THEN
+
+      nbnonnuls = 0
+      do  j =1, coeffraf
+          do i=1, np
+               IF (y(indchildaverage(i,dir) + j -1) .NE. 
+     &               Agrif_SpecialValueFineGrid) THEN
+               nbnonnuls(i) = nbnonnuls(i) + 1
+               x(i) = x(i) +  y(indchildaverage(i,dir) + j -1 )
+               ENDIF
+          End Do
+      enddo
+      do i=1,np
+      x(i)= x(i)*invcoeff(nbnonnuls(i)) 
+      enddo
+
+      ELSE
+!CDIR NOLOOPCHG
+      do  j =1, coeffraf
+!CDIR VECTOR
+          do i=1, np
+             x(i) = x(i) + y(indchildaverage(i,dir) + j -1 )
+          enddo
+      enddo
+       x = x * invcoeffraf
+      ENDIF
+
+      Return 
+C            
+C     
+      End Subroutine average1daftercompute            
 C
 C
