Index: vendors/AGRIF/CMEMS_2020/AGRIF_FILES/modinterpbasic.F90
===================================================================
--- vendors/AGRIF/CMEMS_2020/AGRIF_FILES/modinterpbasic.F90	(revision 10083)
+++ vendors/AGRIF/CMEMS_2020/AGRIF_FILES/modinterpbasic.F90	(revision 10087)
@@ -41,4 +41,5 @@
     integer, dimension(:), allocatable      :: indparentppm_1d, indchildppm_1d
 !
+
     private :: Agrif_limiter_vanleer
 !
@@ -56,12 +57,12 @@
     integer,             intent(in)     :: np           !< Length of input array
     integer,             intent(in)     :: nc           !< Length of output array
-    real,                intent(in)     :: s_parent     !< Parent grid position (s_root = 0)
-    real,                intent(in)     :: s_child      !< Child  grid position (s_root = 0)
-    real,                intent(in)     :: ds_parent    !< Parent grid dx (ds_root = 1)
-    real,                intent(in)     :: ds_child     !< Child  grid dx (ds_root = 1)
+    real(kind=8),                intent(in)     :: s_parent     !< Parent grid position (s_root = 0)
+    real(kind=8),                intent(in)     :: s_child      !< Child  grid position (s_root = 0)
+    real(kind=8),                intent(in)     :: ds_parent    !< Parent grid dx (ds_root = 1)
+    real(kind=8),                intent(in)     :: ds_child     !< Child  grid dx (ds_root = 1)
 !
     integer :: i, coeffraf, locind_parent_left
-    real    :: globind_parent_left, globind_parent_right
-    real    :: invds, invds2, ypos, ypos2, diff
+    real(kind=8)    :: globind_parent_left, globind_parent_right
+    real(kind=8)    :: invds, invds2, ypos, ypos2, diff
 !
     coeffraf = nint(ds_parent/ds_child)
@@ -92,5 +93,9 @@
 !
         diff = globind_parent_right - ypos2
+! quick fix for roundoff error
+        diff=nint(diff*coeffraf)/real(coeffraf)
+
         y(i) = (diff*x(locind_parent_left) + (1.-diff)*x(locind_parent_left+1))
+
         ypos2 = ypos2 + invds2
 !
@@ -104,7 +109,13 @@
     else
         globind_parent_left = s_parent + (locind_parent_left - 1)*ds_parent
-        y(nc) = ((globind_parent_left + ds_parent - ypos)*x(locind_parent_left)  &
-                           + (ypos - globind_parent_left)*x(locind_parent_left+1))*invds
-    endif
+        diff=(globind_parent_left + ds_parent - ypos)*invds
+
+! quick fix for roundoff error
+        diff=nint(diff*coeffraf)/real(coeffraf)
+!        y(nc) = ((globind_parent_left + ds_parent - ypos)*x(locind_parent_left)  &
+!                           + (ypos - globind_parent_left)*x(locind_parent_left+1))*invds
+        y(nc) = (diff*x(locind_parent_left) + (1.-diff)*x(locind_parent_left+1))
+    endif
+
 !---------------------------------------------------------------------------------------------------
 end subroutine Agrif_basicinterp_linear1D
@@ -120,6 +131,6 @@
 !---------------------------------------------------------------------------------------------------
     integer, intent(in) :: np,nc,np2
-    real,    intent(in) :: s_parent, s_child
-    real,    intent(in) :: ds_parent, ds_child
+    real(kind=8),    intent(in) :: s_parent, s_child
+    real(kind=8),    intent(in) :: ds_parent, ds_child
     integer, intent(in) :: dir
 !
@@ -127,7 +138,7 @@
     integer, dimension(:,:), allocatable :: indparent_tmp
     real, dimension(:,:), allocatable :: coeffparent_tmp
-    real    :: ypos,globind_parent_left,globind_parent_right
-    real    :: invds, invds2, invds3
-    real :: ypos2,diff
+    real(kind=8)    :: ypos,globind_parent_left,globind_parent_right
+    real(kind=8)    :: invds, invds2, invds3
+    real(kind=8) :: ypos2,diff
 !
     coeffraf = nint(ds_parent/ds_child)
@@ -164,5 +175,5 @@
         if (ypos2 > globind_parent_right) then
             locind_parent_left = locind_parent_left + 1
-            globind_parent_right = globind_parent_right + 1.
+            globind_parent_right = globind_parent_right + 1.d0
             ypos2 = ypos*invds+(i-1)*invds2
         endif
@@ -239,10 +250,10 @@
     real, dimension(np), intent(in)     :: x
     real, dimension(nc), intent(out)    :: y
-    real,                intent(in)     :: s_parent, s_child
-    real,                intent(in)     :: ds_parent, ds_child
+    real(kind=8),                intent(in)     :: s_parent, s_child
+    real(kind=8),                intent(in)     :: ds_parent, ds_child
 !
     integer :: i, coeffraf, locind_parent_left
-    real    :: ypos,globind_parent_left
-    real    :: deltax, invdsparent
+    real(kind=8)    :: ypos,globind_parent_left
+    real(kind=8)    :: deltax, invdsparent
     real    :: t2,t3,t4,t5,t6,t7,t8
 !
@@ -304,9 +315,9 @@
     real, dimension(np), intent(in)     :: x
     real, dimension(nc), intent(out)    :: y
-    real,                intent(in)     :: s_parent, s_child
-    real,                intent(in)     :: ds_parent, ds_child
+    real(kind=8),                intent(in)     :: s_parent, s_child
+    real(kind=8),                intent(in)     :: ds_parent, ds_child
 !
     integer :: i, coeffraf, locind_parent
-    real    :: ypos
+    real(kind=8)    :: ypos
 !
     coeffraf = nint(ds_parent/ds_child)
@@ -342,10 +353,10 @@
     real, dimension(np), intent(in)     :: x
     real, dimension(nc), intent(out)    :: y
-    real,                intent(in)     :: s_parent, s_child
-    real,                intent(in)     :: ds_parent, ds_child
+    real(kind=8),                intent(in)     :: s_parent, s_child
+    real(kind=8),                intent(in)     :: ds_parent, ds_child
 !
     real, dimension(:), allocatable :: ytemp
     integer :: i,coeffraf,locind_parent_left,locind_parent_last
-    real    :: ypos,xdiffmod,xpmin,xpmax,slope
+    real(kind=8)    :: ypos,xdiffmod,xpmin,xpmax,slope
     integer :: i1,i2,ii
     integer :: diffmod
@@ -429,10 +440,10 @@
     real, dimension(np), intent(in)     :: x
     real, dimension(nc), intent(out)    :: y
-    real,                intent(in)     :: s_parent, s_child
-    real,                intent(in)     :: ds_parent, ds_child
+    real(kind=8),                intent(in)     :: s_parent, s_child
+    real(kind=8),                intent(in)     :: ds_parent, ds_child
 !
     real, dimension(:), allocatable :: ytemp
     integer :: i,coeffraf,locind_parent_left,locind_parent_last
-    real    :: ypos,xdiffmod,xpmin,xpmax,slope
+    real(kind=8)    :: ypos,xdiffmod,xpmin,xpmax,slope
     integer :: i1,i2,ii
     integer :: diffmod
@@ -524,12 +535,13 @@
     real, dimension(np), intent(in)     :: x
     real, dimension(nc), intent(out)    :: y
-    real,                intent(in)     :: s_parent, s_child
-    real,                intent(in)     :: ds_parent, ds_child
+    real(kind=8),                intent(in)     :: s_parent, s_child
+    real(kind=8),                intent(in)     :: ds_parent, ds_child
 !
     integer :: i,coeffraf,locind_parent_left,locind_parent_last
     integer :: iparent,ipos,pos,nmin,nmax
-    real    :: ypos
+    real(kind=8)    :: ypos
     integer :: i1,jj
-    real :: xpmin,a
+    real(kind=8) :: xpmin
+    real :: a
 !
     real, dimension(np) :: xl,delta,a6,slope
@@ -646,6 +658,6 @@
 !---------------------------------------------------------------------------------------------------
     integer,             intent(in)     :: np2, np, nc
-    real,                intent(in)     :: s_parent, s_child
-    real,                intent(in)     :: ds_parent, ds_child
+    real(kind=8),                intent(in)     :: s_parent, s_child
+    real(kind=8),                intent(in)     :: ds_parent, ds_child
     integer,             intent(in)     :: dir
 !
@@ -655,6 +667,7 @@
     integer :: iparent,ipos,pos
     real    :: ypos
-    integer :: i1,jj
-    real :: xpmin,a
+    integer :: i1,jj,k,l,j
+    real(kind=8) :: xpmin
+    real :: a
 !
     integer :: diffmod
@@ -738,7 +751,14 @@
     enddo
 !
-    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)
+
+    k=1
+    l=0
+    do i=1,np2
+      do j=1,nc
+        indchildppm(k,dir) = indchildppm_1d(j)
+        indparentppm(k,dir) = indparentppm_1d(j) + l
+        k=k+1
+      enddo
+      l=l+np
     enddo
 !---------------------------------------------------------------------------------------------------
@@ -746,112 +766,5 @@
 !===================================================================================================
 !
-!===================================================================================================
-!subroutine PPM1dPrecompute(np,nc,&
-!                    s_parent,s_child,ds_parent,ds_child)
-!!
-!!CC   Description:
-!!CC   subroutine to compute coefficient and index for  a 1D interpolation
-!!CC   using piecewise parabolic method
-!!C    Method:
-!!
-!!     Declarations:
-!!
-!      Implicit none
-!!
-!!     Arguments
-!      Integer             :: np,nc
-!!      Real, Dimension(:),Allocatable :: ytemp
-!      Real                :: s_parent,s_child,ds_parent,ds_child
-!!
-!!     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
-!!
-!      Real :: xrmin,xrmax,am3,s2,s1
-!      Real, Dimension(np) :: xl,delta,a6,slope
-!!      Real, Dimension(:),Allocatable  :: diff,diff2,diff3
-!      INTEGER :: diffmod
-!      REAL :: invcoeffraf
-!!
-!      coeffraf = nint(ds_parent/ds_child)
-!!
-!      If (coeffraf == 1) Then
-!          return
-!      End If
-!      invcoeffraf = ds_child/ds_parent
-!!
-!
-!      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
-!!
-!      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)
-!!
-!      xpmin = s_parent + (locind_parent_left-1)*ds_parent
-!      i1 = 1+agrif_int((xpmin-s_child)/ds_child)
-!!
-!!
-!
-!      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
-!!
-!!
-!        diffmod = 0
-!       IF (mod(coeffraf,2) == 0) diffmod = 1
-!!
-!        ipos = i1
-!!
-!        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
-!!
-!        End do
-!
-!      Return
-!      End subroutine ppm1dprecompute
+
 !===================================================================================================
 !
@@ -863,22 +776,32 @@
 ! Use precomputed coefficient and index.
 !---------------------------------------------------------------------------------------------------
-subroutine PPM1dAfterCompute ( x, y, np, nc, dir )
-!---------------------------------------------------------------------------------------------------
+subroutine PPM1dAfterCompute ( x, y, np, nc, dir, indchildppmloc, tabppmloc )
+!---------------------------------------------------------------------------------------------------
+    integer,             intent(in)     :: np, nc
     real, dimension(np), intent(in)     :: x
     real, dimension(nc), intent(out)    :: y
-    integer,             intent(in)     :: np, nc
     integer,             intent(in)     :: dir
-!
-    integer :: i
-!
+    integer, dimension(1:),intent(in) :: indchildppmloc
+    real, dimension(1:,1:),intent(in) :: tabppmloc
+!
+    integer :: i,j,jp1,jp2,jp3,jp4,k
+
+    
     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
+        j = indparentppm(i,dir)
+        jp1=j+1
+        jp2=jp1+1
+        jp3=jp2+1
+        jp4=jp3+1
+        y(i) = tabppmloc(1,indchildppmloc(i)) * x(j) + &
+               tabppmloc(2,indchildppmloc(i)) * x(jp1) + &
+               tabppmloc(3,indchildppmloc(i)) * x(jp2) + &
+               tabppmloc(4,indchildppmloc(i)) * x(jp3) + &
+               tabppmloc(5,indchildppmloc(i)) * x(jp4)
+    enddo
+
 !---------------------------------------------------------------------------------------------------
 end subroutine PPM1dAfterCompute
+
 !===================================================================================================
 !
@@ -1069,13 +992,13 @@
     real, dimension(np), intent(in)     :: x
     real, dimension(nc), intent(out)    :: y
-    real,                intent(in)     :: s_parent, s_child
-    real,                intent(in)     :: ds_parent, ds_child
+    real(kind=8),                intent(in)     :: s_parent, s_child
+    real(kind=8),                intent(in)     :: ds_parent, ds_child
 !
     real, dimension(:), allocatable :: ytemp
     integer :: i,coeffraf,locind_parent_left,locind_parent_last
     integer :: iparent,ipos,pos,nmin,nmax
-    real    :: ypos
+    real(kind=8)    :: ypos
     integer :: i1,jj
-    real :: xpmin
+    real(kind=8) :: xpmin
 !
     real, dimension(np) :: slope
@@ -1166,12 +1089,12 @@
     real, dimension(np), intent(in)     :: x
     real, dimension(nc), intent(out)    :: y
-    real,                intent(in)     :: s_parent, s_child
-    real,                intent(in)     :: ds_parent, ds_child
+    real(kind=8),                intent(in)     :: s_parent, s_child
+    real(kind=8),                intent(in)     :: ds_parent, ds_child
 !
     integer :: i,coeffraf,locind_parent_left,locind_parent_last
     integer :: ipos, pos
-    real    :: ypos,xi
+    real(kind=8)    :: ypos,xi
     integer :: i1,jj
-    real :: xpmin
+    real(kind=8) :: xpmin
 !
     real, dimension(:),   allocatable  :: ytemp
@@ -1276,12 +1199,13 @@
       Real, Dimension(nc) :: y
       Real, Dimension(:),Allocatable :: ytemp
-      Real                :: s_parent,s_child,ds_parent,ds_child
+      Real(kind=8)        :: s_parent,s_child,ds_parent,ds_child
 !
 !     Local scalars
       Integer :: i,coeffraf,locind_parent_left,locind_parent_last
       Integer :: iparent,ipos,pos,nmin,nmax
-      Real    :: ypos
+      Real(kind=8)    :: ypos
       integer :: i1,jj
-      Real :: xpmin,cavg,a,b
+      Real(kind=8) :: xpmin
+      real :: cavg,a,b
 !      
       Real :: xrmin,xrmax,am3,s2,s1  
