Index: branches/2017/dev_r8126_UKMO_AGRIF_vert_interp/NEMOGCM/NEMO/NST_SRC/agrif_oce.F90
===================================================================
--- branches/2017/dev_r8126_UKMO_AGRIF_vert_interp/NEMOGCM/NEMO/NST_SRC/agrif_oce.F90	(revision 8135)
+++ branches/2017/dev_r8126_UKMO_AGRIF_vert_interp/NEMOGCM/NEMO/NST_SRC/agrif_oce.F90	(revision 8855)
@@ -105,221 +105,216 @@
    END FUNCTION agrif_oce_alloc
 
-      subroutine reconstructandremap(tabin,hin,tabout,hout,N,Nout)      
-      implicit none
-      integer N, Nout
-      real tabin(N), tabout(Nout)
-      real hin(N), hout(Nout)
-      real coeffremap(N,3),zwork(N,3)
-      real zwork2(N+1,3)
-      integer k
-      double precision, parameter :: dsmll=1.0d-8  
-      real q,q01,q02,q001,q002,q0
-      real z_win(1:N+1), z_wout(1:Nout+1)
-      real,parameter :: dpthin = 1.D-3
-      integer :: k1, kbox, ktop, ka, kbot
-      real :: tsum, qbot, rpsum, zbox, ztop, zthk, zbot, offset, qtop
+   SUBROUTINE reconstructandremap(tabin,hin,tabout,hout,N,Nout)      
+      !!----------------------------------------------------------------------
+      !!                ***  FUNCTION reconstructandremap  ***
+      !!----------------------------------------------------------------------
+      IMPLICIT NONE
+      INTEGER N, Nout
+      REAL(wp) tabin(N), tabout(Nout)
+      REAL(wp) hin(N), hout(Nout)
+      REAL(wp) coeffremap(N,3),zwork(N,3)
+      REAL(wp) zwork2(N+1,3)
+      INTEGER jk
+      DOUBLE PRECISION, PARAMETER :: dsmll=1.0d-8  
+      REAL(wp) q,q01,q02,q001,q002,q0
+      REAL(wp) z_win(1:N+1), z_wout(1:Nout+1)
+      REAL(wp),PARAMETER :: dpthin = 1.D-3
+      INTEGER :: k1, kbox, ktop, ka, kbot
+      REAL(wp) :: tsum, qbot, rpsum, zbox, ztop, zthk, zbot, offset, qtop
 
       z_win(1)=0.; z_wout(1)= 0.
-      do k=1,N
-       z_win(k+1)=z_win(k)+hin(k)
-      enddo 
+      DO jk=1,N
+         z_win(jk+1)=z_win(jk)+hin(jk)
+      ENDDO 
       
-      do k=1,Nout
-       z_wout(k+1)=z_wout(k)+hout(k)       
-      enddo       
-
-        do k=2,N
-          zwork(k,1)=1./(hin(k-1)+hin(k))
-        enddo
-        
-        do k=2,N-1
-          q0 = 1./(hin(k-1)+hin(k)+hin(k+1))
-          zwork(k,2)=hin(k-1)+2.*hin(k)+hin(k+1)
-          zwork(k,3)=q0
-        enddo       
-      
-        do k= 2,N
-        zwork2(k,1)=zwork(k,1)*(tabin(k)-tabin(k-1))
-        enddo
-        
-        coeffremap(:,1) = tabin(:)
+      DO jk=1,Nout
+         z_wout(jk+1)=z_wout(jk)+hout(jk)       
+      ENDDO       
+
+      DO jk=2,N
+         zwork(jk,1)=1./(hin(jk-1)+hin(jk))
+      ENDDO
+        
+      DO jk=2,N-1
+         q0 = 1./(hin(jk-1)+hin(jk)+hin(jk+1))
+         zwork(jk,2)=hin(jk-1)+2.*hin(jk)+hin(jk+1)
+         zwork(jk,3)=q0
+      ENDDO       
+     
+      DO jk= 2,N
+         zwork2(jk,1)=zwork(jk,1)*(tabin(jk)-tabin(jk-1))
+      ENDDO
+        
+      coeffremap(:,1) = tabin(:)
  
-         do k=2,N-1
-        q001 = hin(k)*zwork2(k+1,1)
-        q002 = hin(k)*zwork2(k,1)        
-        if (q001*q002 < 0) then
-          q001 = 0.
-          q002 = 0.
-        endif
-        q=zwork(k,2)
-        q01=q*zwork2(k+1,1)
-        q02=q*zwork2(k,1)
-        if (abs(q001) > abs(q02)) q001 = q02
-        if (abs(q002) > abs(q01)) q002 = q01
-        
-        q=(q001-q002)*zwork(k,3)
-        q001=q001-q*hin(k+1)
-        q002=q002+q*hin(k-1)
-        
-        coeffremap(k,3)=coeffremap(k,1)+q001
-        coeffremap(k,2)=coeffremap(k,1)-q002
-        
-        zwork2(k,1)=(2.*q001-q002)**2
-        zwork2(k,2)=(2.*q002-q001)**2
-        enddo
-        
-        do k=1,N
-        if     (k.eq.1 .or. k.eq.N .or.   hin(k).le.dpthin) then
-        coeffremap(k,3) = coeffremap(k,1)
-        coeffremap(k,2) = coeffremap(k,1)
-        zwork2(k,1) = 0.
-        zwork2(k,2) = 0.
-        endif
-        enddo
-        
-        do k=2,N
-        q002=max(zwork2(k-1,2),dsmll)
-        q001=max(zwork2(k,1),dsmll)
-        zwork2(k,3)=(q001*coeffremap(k-1,3)+q002*coeffremap(k,2))/(q001+q002)
-        enddo
-        
-        zwork2(1,3) = 2*coeffremap(1,1)-zwork2(2,3)
-        zwork2(N+1,3)=2*coeffremap(N,1)-zwork2(N,3)
+      DO jk=2,N-1
+         q001 = hin(jk)*zwork2(jk+1,1)
+         q002 = hin(jk)*zwork2(jk,1)        
+         IF (q001*q002 < 0) then
+            q001 = 0.
+            q002 = 0.
+         ENDIF
+         q=zwork(jk,2)
+         q01=q*zwork2(jk+1,1)
+         q02=q*zwork2(jk,1)
+         IF (abs(q001) > abs(q02)) q001 = q02
+         IF (abs(q002) > abs(q01)) q002 = q01
+        
+         q=(q001-q002)*zwork(jk,3)
+         q001=q001-q*hin(jk+1)
+         q002=q002+q*hin(jk-1)
+        
+         coeffremap(jk,3)=coeffremap(jk,1)+q001
+         coeffremap(jk,2)=coeffremap(jk,1)-q002
+        
+         zwork2(jk,1)=(2.*q001-q002)**2
+         zwork2(jk,2)=(2.*q002-q001)**2
+      ENDDO
+        
+      DO jk=1,N
+         IF(jk.EQ.1 .OR. jk.EQ.N .OR. hin(jk).LE.dpthin) THEN
+            coeffremap(jk,3) = coeffremap(jk,1)
+            coeffremap(jk,2) = coeffremap(jk,1)
+            zwork2(jk,1) = 0.
+            zwork2(jk,2) = 0.
+         ENDIF
+      ENDDO
+        
+      DO jk=2,N
+         q002=max(zwork2(jk-1,2),dsmll)
+         q001=max(zwork2(jk,1),dsmll)
+         zwork2(jk,3)=(q001*coeffremap(jk-1,3)+q002*coeffremap(jk,2))/(q001+q002)
+      ENDDO
+        
+      zwork2(1,3) = 2*coeffremap(1,1)-zwork2(2,3)
+      zwork2(N+1,3)=2*coeffremap(N,1)-zwork2(N,3)
  
-        do k=1,N
-        q01=zwork2(k+1,3)-coeffremap(k,1)
-        q02=coeffremap(k,1)-zwork2(k,3)
-        q001=2.*q01
-        q002=2.*q02
-        if (q01*q02<0) then
-          q01=0.
-          q02=0.
-        elseif (abs(q01)>abs(q002)) then
-          q01=q002
-        elseif (abs(q02)>abs(q001)) then
-          q02=q001
-        endif
-        coeffremap(k,2)=coeffremap(k,1)-q02
-        coeffremap(k,3)=coeffremap(k,1)+q01
-        enddo
+      DO jk=1,N
+         q01=zwork2(jk+1,3)-coeffremap(jk,1)
+         q02=coeffremap(jk,1)-zwork2(jk,3)
+         q001=2.*q01
+         q002=2.*q02
+         IF (q01*q02<0) then
+            q01=0.
+            q02=0.
+         ELSEIF (abs(q01)>abs(q002)) then
+            q01=q002
+         ELSEIF (abs(q02)>abs(q001)) then
+            q02=q001
+         ENDIF
+         coeffremap(jk,2)=coeffremap(jk,1)-q02
+         coeffremap(jk,3)=coeffremap(jk,1)+q01
+      ENDDO
 
       zbot=0.0
       kbot=1
-      do k=1,Nout
-        ztop=zbot  !top is bottom of previous layer
-        ktop=kbot
-        if     (ztop.ge.z_win(ktop+1)) then
-          ktop=ktop+1
-        endif
-        
-        zbot=z_wout(k+1)
-        zthk=zbot-ztop
-
-        if     (zthk.gt.dpthin .and. ztop.lt.z_wout(Nout+1)) then
-
-          kbot=ktop
-          do while (z_win(kbot+1).lt.zbot.and.kbot.lt.N)
-            kbot=kbot+1
-          enddo
-          zbox=zbot
-          do k1= k+1,Nout
-            if     (z_wout(k1+1)-z_wout(k1).gt.dpthin) then
-              exit !thick layer
-            else
-              zbox=z_wout(k1+1)  !include thin adjacent layers
-              if     (zbox.eq.z_wout(Nout+1)) then
-                exit !at bottom
-              endif
-            endif
-          enddo
-          zthk=zbox-ztop
-
-          kbox=ktop
-          do while (z_win(kbox+1).lt.zbox.and.kbox.lt.N)
-            kbox=kbox+1
-          enddo
+      DO jk=1,Nout
+         ztop=zbot  !top is bottom of previous layer
+         ktop=kbot
+         IF     (ztop.GE.z_win(ktop+1)) then
+            ktop=ktop+1
+         ENDIF
+        
+         zbot=z_wout(jk+1)
+         zthk=zbot-ztop
+
+         IF(zthk.GT.dpthin .AND. ztop.LT.z_wout(Nout+1)) THEN
+
+            kbot=ktop
+            DO while (z_win(kbot+1).lt.zbot.and.kbot.lt.N)
+               kbot=kbot+1
+            ENDDO
+            zbox=zbot
+            DO k1= jk+1,Nout
+               IF     (z_wout(k1+1)-z_wout(k1).GT.dpthin) THEN
+                  exit !thick layer
+               ELSE
+                  zbox=z_wout(k1+1)  !include thin adjacent layers
+                  IF(zbox.EQ.z_wout(Nout+1)) THEN
+                     exit !at bottom
+                  ENDIF
+               ENDIF
+            ENDDO
+            zthk=zbox-ztop
+
+            kbox=ktop
+            DO while (z_win(kbox+1).lt.zbox.and.kbox.lt.N)
+               kbox=kbox+1
+            ENDDO
           
-          if     (ktop.eq.kbox) then
-
-
-            if     (z_wout(k)  .ne.z_win(kbox)   .or.z_wout(k+1).ne.z_win(kbox+1)     ) then
-
-              if     (hin(kbox).gt.dpthin) then
-                q001 = (zbox-z_win(kbox))/hin(kbox)
-                q002 = (ztop-z_win(kbox))/hin(kbox)
-                q01=q001**2+q002**2+q001*q002+1.-2.*(q001+q002)
-                q02=q01-1.+(q001+q002)
-                q0=1.-q01-q02
-              else
-                q0 = 1.0
-                q01 = 0.
-                q02 = 0.
-              endif
-          tabout(k)=q0*coeffremap(kbox,1)+q01*coeffremap(kbox,2)+q02*coeffremap(kbox,3)
+            IF(ktop.EQ.kbox) THEN
+               IF(z_wout(jk).NE.z_win(kbox).OR.z_wout(jk+1).NE.z_win(kbox+1)) THEN
+                  IF(hin(kbox).GT.dpthin) THEN
+                     q001 = (zbox-z_win(kbox))/hin(kbox)
+                     q002 = (ztop-z_win(kbox))/hin(kbox)
+                     q01=q001**2+q002**2+q001*q002+1.-2.*(q001+q002)
+                     q02=q01-1.+(q001+q002)
+                     q0=1.-q01-q02
+                  ELSE
+                     q0 = 1.0
+                     q01 = 0.
+                     q02 = 0.
+                  ENDIF
+                  tabout(jk)=q0*coeffremap(kbox,1)+q01*coeffremap(kbox,2)+q02*coeffremap(kbox,3)
+               ELSE
+                  tabout(jk) = tabin(kbox)
+               ENDIF 
+            ELSE
+               IF(ktop.LE.jk .AND. kbox.GE.jk) THEN
+                  ka = jk
+               ELSEIF (kbox-ktop.GE.3) THEN
+                  ka = (kbox+ktop)/2
+               ELSEIF (hin(ktop).GE.hin(kbox)) THEN
+                  ka = ktop
+               ELSE
+                  ka = kbox
+               ENDIF !choose ka
+    
+               offset=coeffremap(ka,1)
+    
+               qtop = z_win(ktop+1)-ztop !partial layer thickness
+               IF(hin(ktop).GT.dpthin) THEN
+                  q=(ztop-z_win(ktop))/hin(ktop)
+                  q01=q*(q-1.)
+                  q02=q01+q
+                  q0=1-q01-q02            
+               ELSE
+                  q0 = 1.
+                  q01 = 0.
+                  q02 = 0.
+               ENDIF
+               
+               tsum =((q0*coeffremap(ktop,1)+q01*coeffremap(ktop,2)+q02*coeffremap(ktop,3))-offset)*qtop
+    
+               DO k1= ktop+1,kbox-1
+                  tsum =tsum +(coeffremap(k1,1)-offset)*hin(k1)
+               ENDDO !k1
+               
+               qbot = zbox-z_win(kbox) !partial layer thickness
+               IF(hin(kbox).GT.dpthin) THEN
+                  q=qbot/hin(kbox)
+                  q01=(q-1.)**2
+                  q02=q01-1.+q
+                  q0=1-q01-q02                            
+               ELSE
+                  q0 = 1.0
+                  q01 = 0.
+                  q02 = 0.
+               ENDIF
               
-            else
-            tabout(k) = tabin(kbox)
-              
-            endif 
-
-          else
-
-            if     (ktop.le.k .and. kbox.ge.k) then
-              ka = k
-            elseif (kbox-ktop.ge.3) then
-              ka = (kbox+ktop)/2
-            elseif (hin(ktop).ge.hin(kbox)) then
-              ka = ktop
-            else
-              ka = kbox
-            endif !choose ka
-
-            offset=coeffremap(ka,1)
-
-            qtop = z_win(ktop+1)-ztop !partial layer thickness
-            if     (hin(ktop).gt.dpthin) then
-              q=(ztop-z_win(ktop))/hin(ktop)
-              q01=q*(q-1.)
-              q02=q01+q
-              q0=1-q01-q02            
-            else
-              q0 = 1.
-              q01 = 0.
-              q02 = 0.
-            endif
-            
-            tsum =((q0*coeffremap(ktop,1)+q01*coeffremap(ktop,2)+q02*coeffremap(ktop,3))-offset)*qtop
-
-            do k1= ktop+1,kbox-1
-              tsum =tsum +(coeffremap(k1,1)-offset)*hin(k1)
-            enddo !k1
-
-            
-            qbot = zbox-z_win(kbox) !partial layer thickness
-            if     (hin(kbox).gt.dpthin) then
-              q=qbot/hin(kbox)
-              q01=(q-1.)**2
-              q02=q01-1.+q
-              q0=1-q01-q02                            
-            else
-              q0 = 1.0
-              q01 = 0.
-              q02 = 0.
-            endif
-           
-            tsum = tsum +((q0*coeffremap(kbox,1)+q01*coeffremap(kbox,2)+q02*coeffremap(kbox,3))-offset)*qbot
-            
-            rpsum=1.0d0/zthk
-              tabout(k)=offset+tsum*rpsum
-              
-          endif !single or multiple layers
-        else
-        if (k==1) then
-        write(*,'(a7,i4,i4,3f12.5)')'problem = ',N,Nout,zthk,z_wout(k+1),hout(1)
-        endif
-         tabout(k) = tabout(k-1)
-          
-        endif !normal:thin layer
-      enddo !k
+               tsum = tsum +((q0*coeffremap(kbox,1)+q01*coeffremap(kbox,2)+q02*coeffremap(kbox,3))-offset)*qbot
+               
+               rpsum=1.0d0/zthk
+               tabout(jk)=offset+tsum*rpsum
+                 
+            ENDIF !single or multiple layers
+         ELSE
+            IF (jk==1) THEN
+               write(*,'(a7,i4,i4,3f12.5)')'problem = ',N,Nout,zthk,z_wout(jk+1),hout(1)
+            ENDIF
+            tabout(jk) = tabout(jk-1)
+             
+         ENDIF !normal:thin layer
+      ENDDO !jk
             
       return
