Index: /NEMO/branches/2020/dev_r13312_AGRIF-03-04_jchanut_vinterp_tstep/src/NST/agrif_top_update.F90
===================================================================
--- /NEMO/branches/2020/dev_r13312_AGRIF-03-04_jchanut_vinterp_tstep/src/NST/agrif_top_update.F90	(revision 14047)
+++ /NEMO/branches/2020/dev_r13312_AGRIF-03-04_jchanut_vinterp_tstep/src/NST/agrif_top_update.F90	(revision 14048)
@@ -40,6 +40,8 @@
       IF (Agrif_Root()) RETURN 
       !
-      Agrif_UseSpecialValueInUpdate = .TRUE.
+      l_vremap                      = ln_vert_remap
+      Agrif_UseSpecialValueInUpdate = .NOT.l_vremap
       Agrif_SpecialValueFineGrid    = 0._wp
+
       ! 
 # if ! defined DECAL_FEEDBACK
@@ -52,12 +54,10 @@
       !
       Agrif_UseSpecialValueInUpdate = .FALSE.
+      l_vremap                      = .FALSE.
       !
    END SUBROUTINE Agrif_Update_Trc
 
-#ifdef key_vertical
    SUBROUTINE updateTRC( tabres, i1, i2, j1, j2, k1, k2, n1, n2, before )
-      !!---------------------------------------------
-      !!           *** ROUTINE updateT ***
-      !!---------------------------------------------
+
       INTEGER, INTENT(in) :: i1,i2,j1,j2,k1,k2,n1,n2
       REAL(wp),DIMENSION(i1:i2,j1:j2,k1:k2,n1:n2), INTENT(inout) :: tabres
@@ -72,169 +72,131 @@
       REAL(wp) :: tabin(k1:k2,1:jptra)
       REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk,1:jptra) :: tabres_child
-      !!---------------------------------------------
-      !
+
       IF (before) THEN
-         AGRIF_SpecialValue = -999._wp
-         DO jn = n1,n2-1
+         IF ( l_vremap ) THEN
+            DO jn = n1,n2-1
+               DO jk=k1,k2
+                  DO jj=j1,j2
+                     DO ji=i1,i2
+                        tabres(ji,jj,jk,jn) = tr(ji,jj,jk,jn,Kmm_a) * e3t(ji,jj,jk,Kmm_a)
+                     END DO
+                  END DO
+               END DO
+            END DO
             DO jk=k1,k2
                DO jj=j1,j2
                   DO ji=i1,i2
-                     tabres(ji,jj,jk,jn) = (tr(ji,jj,jk,jn,Kmm_a) * e3t(ji,jj,jk,Kmm_a) ) &
-                                           * tmask(ji,jj,jk) + (tmask(ji,jj,jk)-1)*999._wp
-                  END DO
-               END DO
-            END DO
-         END DO
-         DO jk=k1,k2
+                     tabres(ji,jj,jk,n2) = tmask(ji,jj,jk) * e3t(ji,jj,jk,Kmm_a)
+                  END DO
+               END DO
+            END DO
+         ELSE
+            DO jn = 1,jptra
+               DO jk=k1,k2
+                  DO jj=j1,j2
+                     DO ji=i1,i2
+                        tabres(ji,jj,jk,jn) = tr(ji,jj,jk,jn,Kmm_a)  * e3t(ji,jj,jk,Kmm_a) / e3t_0(ji,jj,jk)
+                     END DO
+                  END DO
+               END DO
+            END DO
+
+         ENDIF
+      ELSE
+         IF ( l_vremap ) THEN
+            tabres_child(:,:,:,:) = 0._wp
+            AGRIF_SpecialValue = 0._wp
             DO jj=j1,j2
                DO ji=i1,i2
-                  tabres(ji,jj,jk,n2) = tmask(ji,jj,jk) * e3t(ji,jj,jk,Kmm_a) &
-                                           + (tmask(ji,jj,jk)-1)*999._wp
-               END DO
-            END DO
-         END DO
-      ELSE
-         tabres_child(:,:,:,:) = 0.
-         AGRIF_SpecialValue = 0._wp
-         DO jj=j1,j2
-            DO ji=i1,i2
-               N_in = 0
-               DO jk=k1,k2 !k2 = jpk of child grid
-                  IF (tabres(ji,jj,jk,n2) == 0  ) EXIT
-                  N_in = N_in + 1
-                  tabin(jk,:) = tabres(ji,jj,jk,n1:n2-1)/tabres(ji,jj,jk,n2)
-                  h_in(N_in) = tabres(ji,jj,jk,n2)
+                  N_in = 0
+                  DO jk=k1,k2 !k2 = jpk of child grid
+                     IF (tabres(ji,jj,jk,n2) <= 1.e-6_wp  ) EXIT
+                     N_in = N_in + 1
+                     tabin(jk,:) = tabres(ji,jj,jk,n1:n2-1)/tabres(ji,jj,jk,n2)
+                     h_in(N_in) = tabres(ji,jj,jk,n2)
+                  ENDDO
+                  N_out = 0
+                  DO jk=1,jpk ! jpk of parent grid
+                     IF (tmask(ji,jj,jk) == 0 ) EXIT ! TODO: Will not work with ISF
+                     N_out = N_out + 1
+                     h_out(N_out) = e3t(ji,jj,jk,Kmm_a)
+                  ENDDO
+                  IF (N_in*N_out > 0) THEN !Remove this?
+                     CALL reconstructandremap(tabin(1:N_in,1:jptra),h_in(1:N_in),tabres_child(ji,jj,1:N_out,1:jptra),h_out(1:N_out),N_in,N_out,jptra)
+                  ENDIF
                ENDDO
-               N_out = 0
-               DO jk=1,jpk ! jpk of parent grid
-                  IF (tmask(ji,jj,jk) < -900) EXIT ! TODO: Will not work with ISF
-                  N_out = N_out + 1
-                  h_out(N_out) = e3t(ji,jj,jk,Kmm_a) !Parent grid scale factors. Could multiply by e1e2t here instead of division above
-               ENDDO
-               IF (N_in > 0) THEN !Remove this?
-                  h_diff = sum(h_out(1:N_out))-sum(h_in(1:N_in))
-                  IF (h_diff < -1.e-4) THEN
-                     print *,'CHECK YOUR bathy T points ...',ji,jj,h_diff,sum(h_in(1:N_in)),sum(h_out(1:N_out))
-                     print *,h_in(1:N_in)
-                     print *,h_out(1:N_out)
-                     STOP
-                  ENDIF
-                  CALL reconstructandremap(tabin(1:N_in,1:jptra),h_in(1:N_in),tabres_child(ji,jj,1:N_out,1:jptra),h_out(1:N_out),N_in,N_out,jptra)
-               ENDIF
             ENDDO
-         ENDDO
-         !
-         IF (.NOT.(lk_agrif_fstep.AND.(l_1st_euler))) THEN
-            ! Add asselin part
-            DO jn = 1,jptra
-               DO jk=1,jpkm1
+
+            IF (.NOT.(lk_agrif_fstep.AND.(l_1st_euler))) THEN
+               ! Add asselin part
+               DO jn = 1,jptra
+                  DO jk = 1, jpkm1
+                     DO jj = j1, j2
+                        DO ji = i1, i2
+                           IF( tabres_child(ji,jj,jk,jn) /= 0._wp ) THEN
+                              ztb  = tr(ji,jj,jk,jn,Kbb_a) * e3t(ji,jj,jk,Kbb_a) ! fse3t_b prior update should be used
+                              ztnu = tabres_child(ji,jj,jk,jn) * e3t(ji,jj,jk,Kmm_a)
+                              ztno = tr(ji,jj,jk,jn,Kmm_a) * e3t(ji,jj,jk,Krhs_a)
+                              tr(ji,jj,jk,jn,Kbb_a) = ( ztb + rn_atfp * ( ztnu - ztno) )  &
+                                        &        * tmask(ji,jj,jk) / e3t(ji,jj,jk,Kbb_a)
+                           ENDIF
+                        END DO
+                     END DO
+                  END DO
+               END DO
+            ENDIF
+            DO jn = 1,jptra
+               DO jk = 1, jpkm1
+                  DO jj = j1, j2
+                     DO ji = i1, i2
+                        IF( tabres_child(ji,jj,jk,jn) /= 0._wp ) THEN
+                           tr(ji,jj,jk,jn,Kmm_a) = tabres_child(ji,jj,jk,jn)
+                        END IF
+                     END DO
+                  END DO
+               END DO
+            END DO
+         ELSE
+            DO jn = 1,jptra
+               tabres(i1:i2,j1:j2,k1:k2,jn) =  tabres(i1:i2,j1:j2,k1:k2,jn) * e3t_0(i1:i2,j1:j2,k1:k2) &
+                                            & * tmask(i1:i2,j1:j2,k1:k2)
+            ENDDO
+            IF (.NOT.(lk_agrif_fstep.AND.(l_1st_euler))) THEN
+               ! Add asselin part
+               DO jn = 1,jptra
+                  DO jk = k1, k2
+                     DO jj = j1, j2
+                        DO ji = i1, i2
+                           IF( tabres(ji,jj,jk,jn) /= 0._wp ) THEN
+                              ztb  = tr(ji,jj,jk,jn,Kbb_a) * e3t(ji,jj,jk,Kbb_a) ! fse3t_b prior update should be used
+                              ztnu = tabres(ji,jj,jk,jn)
+                              ztno = tr(ji,jj,jk,jn,Kmm_a) * e3t(ji,jj,jk,Krhs_a)
+                              tr(ji,jj,jk,jn,Kbb_a) = ( ztb + rn_atfp * ( ztnu - ztno) )  &
+                                        &        * tmask(ji,jj,jk) / e3t(ji,jj,jk,Kbb_a)
+                           ENDIF
+                        END DO
+                     END DO
+                  END DO
+               END DO
+            ENDIF
+            DO jn = 1,jptra
+               DO jk=k1,k2
                   DO jj=j1,j2
                      DO ji=i1,i2
-                        IF( tabres_child(ji,jj,jk,jn) .NE. 0. ) THEN
-                           ztb  = tr(ji,jj,jk,jn,Kbb_a) * e3t(ji,jj,jk,Kbb_a) ! fse3t_b prior update should be used
-                           ztnu = tabres_child(ji,jj,jk,jn) * e3t(ji,jj,jk,Kmm_a)
-                           ztno = tr(ji,jj,jk,jn,Kmm_a) * e3t(ji,jj,jk,Krhs_a)
-                           tr(ji,jj,jk,jn,Kbb_a) = ( ztb + rn_atfp * ( ztnu - ztno) )  & 
-                                     &        * tmask(ji,jj,jk) / e3t(ji,jj,jk,Kbb_a)
-                        ENDIF
-                     ENDDO
-                  ENDDO
-               ENDDO
-            ENDDO
+                        IF( tabres(ji,jj,jk,jn) /= 0._wp ) THEN
+                           tr(ji,jj,jk,jn,Kmm_a) = tabres(ji,jj,jk,jn) / e3t(ji,jj,jk,Kmm_a)
+                        END IF
+                     END DO
+                  END DO
+               END DO
+            END DO
+            !
          ENDIF
-         DO jn = 1,jptra
-            DO jk=1,jpkm1
-               DO jj=j1,j2
-                  DO ji=i1,i2
-                     IF( tabres_child(ji,jj,jk,jn) .NE. 0. ) THEN 
-                        tr(ji,jj,jk,jn,Kmm_a) = tabres_child(ji,jj,jk,jn)
-                     END IF
-                  END DO
-               END DO
-            END DO
-         END DO
-         !
          IF  ((l_1st_euler).AND.(Agrif_Nb_Step()==0) ) THEN
             tr(i1:i2,j1:j2,1:jpkm1,1:jptra,Kbb_a)  = tr(i1:i2,j1:j2,1:jpkm1,1:jptra,Kmm_a)
          ENDIF
-         !
-
       ENDIF
-      ! 
+      !
    END SUBROUTINE updateTRC
-
-
-#else
-   SUBROUTINE updateTRC( tabres, i1, i2, j1, j2, k1, k2, n1, n2, before )
-      !!----------------------------------------------------------------------
-      !!                      *** ROUTINE updateTRC ***
-      !!----------------------------------------------------------------------
-      INTEGER                                    , INTENT(in   ) ::   i1, i2, j1, j2, k1, k2, n1, n2
-      REAL(wp),DIMENSION(i1:i2,j1:j2,k1:k2,n1:n2), INTENT(inout) ::   tabres
-      LOGICAL                                    , INTENT(in   ) ::   before
-      !!
-      INTEGER :: ji,jj,jk,jn
-      REAL(wp) :: ztb, ztnu, ztno
-      !!----------------------------------------------------------------------
-      !
-      !
-      IF (before) THEN
-         DO jn = n1,n2
-            DO jk=k1,k2
-               DO jj=j1,j2
-                  DO ji=i1,i2
-!> jc tmp
-                     tabres(ji,jj,jk,jn) = tr(ji,jj,jk,jn,Kmm_a)  * e3t(ji,jj,jk,Kmm_a) / e3t_0(ji,jj,jk)
-!                     tabres(ji,jj,jk,jn) = tr(ji,jj,jk,jn,Kmm_a)  * e3t(ji,jj,jk,Kmm_a)
-!< jc tmp
-                  END DO
-               END DO
-            END DO
-         END DO
-      ELSE
-!> jc tmp
-         DO jn = n1,n2
-            tabres(i1:i2,j1:j2,k1:k2,jn) =  tabres(i1:i2,j1:j2,k1:k2,jn) * e3t_0(i1:i2,j1:j2,k1:k2) &
-                                         & * tmask(i1:i2,j1:j2,k1:k2)
-         ENDDO
-!< jc tmp
-         IF (.NOT.(lk_agrif_fstep.AND.(l_1st_euler))) THEN
-            ! Add asselin part
-            DO jn = n1,n2
-               DO jk=k1,k2
-                  DO jj=j1,j2
-                     DO ji=i1,i2
-                        IF( tabres(ji,jj,jk,jn) .NE. 0. ) THEN
-                           ztb  = tr(ji,jj,jk,jn,Kbb_a) * e3t(ji,jj,jk,Kbb_a) ! fse3t_b prior update should be used
-                           ztnu = tabres(ji,jj,jk,jn)
-                           ztno = tr(ji,jj,jk,jn,Kmm_a) * e3t(ji,jj,jk,Krhs_a)
-                           tr(ji,jj,jk,jn,Kbb_a) = ( ztb + rn_atfp * ( ztnu - ztno) )  & 
-                                     &        * tmask(ji,jj,jk) / e3t(ji,jj,jk,Kbb_a)
-                        ENDIF
-                     ENDDO
-                  ENDDO
-               ENDDO
-            ENDDO
-         ENDIF
-         DO jn = n1,n2
-            DO jk=k1,k2
-               DO jj=j1,j2
-                  DO ji=i1,i2
-                     IF( tabres(ji,jj,jk,jn) .NE. 0. ) THEN 
-                        tr(ji,jj,jk,jn,Kmm_a) = tabres(ji,jj,jk,jn) / e3t(ji,jj,jk,Kmm_a)
-                     END IF
-                  END DO
-               END DO
-            END DO
-         END DO
-         !
-         IF  ((l_1st_euler).AND.(Agrif_Nb_Step()==0) ) THEN
-            tr(i1:i2,j1:j2,k1:k2,n1:n2,Kbb_a)  = tr(i1:i2,j1:j2,k1:k2,n1:n2,Kmm_a)
-         ENDIF
-         !
-      ENDIF
-      ! 
-   END SUBROUTINE updateTRC
-#endif
 
 #else
