Index: branches/UKMO/AMM15_v3_6_STABLE_package_collate_utils366_fabmv1_v2/NEMOGCM/NEMO/TOP_SRC/FABM/vertical_movement_fabm.F90
===================================================================
--- branches/UKMO/AMM15_v3_6_STABLE_package_collate_utils366_fabmv1_v2/NEMOGCM/NEMO/TOP_SRC/FABM/vertical_movement_fabm.F90	(revision 13488)
+++ branches/UKMO/AMM15_v3_6_STABLE_package_collate_utils366_fabmv1_v2/NEMOGCM/NEMO/TOP_SRC/FABM/vertical_movement_fabm.F90	(revision 13489)
@@ -15,5 +15,4 @@
    USE trc
    USE par_fabm
-   USE fabm
    USE dom_oce
 #if defined key_trdtrc && defined key_iomput
@@ -25,4 +24,5 @@
 
 #  include "domzgr_substitute.h90"
+#  include "vectopt_loop_substitute.h90"
 
    PRIVATE
@@ -32,7 +32,4 @@
    ! Work arrays for vertical advection (residual movement/sinking/floating)
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, TARGET, DIMENSION(:,:,:) :: w_ct
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, TARGET, DIMENSION(:,:)   :: w_if
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, TARGET, DIMENSION(:,:)   :: zwgt_if
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, TARGET, DIMENSION(:,:)   :: flux_if
 #if defined key_trdtrc && defined key_iomput
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, TARGET, DIMENSION(:,:,:,:) :: tr_vmv
@@ -41,23 +38,20 @@
    CONTAINS
 
-   SUBROUTINE compute_vertical_movement( kt )
+   SUBROUTINE compute_vertical_movement( kt, method )
       !!----------------------------------------------------------------------
       !!                     ***  compute_vertical_movement  ***
       !!
-      !! ** Purpose :   compute vertical movement of FABM tracers
+      !! ** Purpose : compute vertical movement of FABM tracers through the water
+      !!              (sinking/floating/active movement)
       !!
-      !! ** Method  : Sets additional vertical velocity field and computes
-      !!              resulting advection using a conservative 3rd upwind
-      !!              scheme with QUICKEST TVD limiter, based on GOTM
-      !!              module adv_center.F90 (www.gotm.net). Currently assuming
-      !!              zero flux at sea surface and sea floor.
+      !! ** Method  : Retrieves additional vertical velocity field and applies
+      !!              advection scheme.
       !!----------------------------------------------------------------------
       !
-      INTEGER, INTENT(in) ::   kt   ! ocean time-step index
-      INTEGER :: ji,jj,jk,jn,k_floor,n_iter,n_count
-      INTEGER,PARAMETER :: n_itermax=100
-      REAL(wp) :: cmax_no,z2dt
-      REAL(wp),DIMENSION(jpk) :: tr_it,tr_u,tr_d,tr_c,tr_slope,c_no,flux_lim
-      REAL(wp),DIMENSION(jpk) :: phi_lim,x_fac
+      INTEGER, INTENT(in) ::   kt     ! ocean time-step index
+      INTEGER, INTENT(in) ::   method ! advection method (1: 1st order upstream, 3: 3rd order TVD with QUICKEST limiter)
+
+      INTEGER :: ji,jj,jk,jn,k_floor
+      REAL(wp) :: zwgt_if(1:jpkm1-1), dc(1:jpkm1), w_if(1:jpkm1-1), z2dt, h(1:jpkm1)
 #if defined key_trdtrc
       CHARACTER (len=20) :: cltra
@@ -69,144 +63,53 @@
 
       IF( neuler == 0 .AND. kt == nittrc000 ) THEN
-          z2dt = rdt                  ! set time step size (Euler)
+         z2dt = rdt                  ! set time step size (Euler)
       ELSE
-          z2dt = 2._wp * rdt          ! set time step size (Leapfrog)
+         z2dt = 2._wp * rdt          ! set time step size (Leapfrog)
       ENDIF
+
       ! Compute interior vertical velocities and include them in source array.
-      DO jj=1,jpj ! j-loop
-         ! Get vertical velocities at layer centres (entire 1:jpi,1:jpk slice).
-         DO jk=1,jpk
-            CALL fabm_get_vertical_movement(model,1,jpi,jj,jk,w_ct(:,jk,:))
+      DO jj=2,jpjm1 ! j-loop
+         ! Get vertical velocities at layer centres (entire i-k slice).
+         DO jk=1,jpkm1
+            CALL model%get_vertical_movement(fs_2,fs_jpim1,jj,jk,w_ct(:,jk,:))
          END DO
-
-         DO ji=1,jpi ! i-loop
+         DO ji=fs_2,fs_jpim1 ! i-loop
             ! Only process this horizontal point (ji,jj) if number of layers exceeds 1
-            IF (mbkt(ji,jj)>1) THEN ! Level check
-               k_floor=mbkt(ji,jj)
+            k_floor = mbkt(ji,jj)
+            IF (k_floor > 1) THEN
                ! Linearly interpolate to velocities at the interfaces between layers
                ! Note:
-               !    - interface k sits between cell centre k and k-1,
-               !    - k [1,jpk] increases downwards
+               !    - interface k sits between cell centre k and k+1 (k=0 for surface)
+               !    - k [1,jpkm1] increases downwards
                !    - upward velocity is positive, downward velocity is negative
-               zwgt_if(1,:)=0._wp ! surface
-               w_if(1,:)=0._wp ! surface
-               zwgt_if(2:k_floor,:)=spread(&
-                   fse3t(ji,jj,2:k_floor)/ (fse3t(ji,jj,1:k_floor-1)+fse3t(ji,jj,2:k_floor))&
-                   ,2,jp_fabm)
-               w_if(2:k_floor,:) = zwgt_if(2:k_floor,:)*w_ct(ji,1:k_floor-1,:)&
-                  +(1._wp-zwgt_if(1:k_floor-1,:))*w_ct(ji,2:k_floor,:)
-               zwgt_if(k_floor+1:,:)=0._wp ! sea floor and below
-               w_if(k_floor+1:,:)=0._wp ! sea floor and below
+               h(1:k_floor) = fse3t(ji,jj,1:k_floor)
+               zwgt_if(1:k_floor-1) = h(2:k_floor) / (h(1:k_floor-1) + h(2:k_floor))
 
                ! Advect:
                DO jn=1,jp_fabm ! State loop
-                  ! get maximum Courant number:
-                  c_no(2:k_floor)=abs(w_if(2:k_floor,jn))*z2dt/ &
-                                ( 0.5_wp*(fse3t(ji,jj,2:k_floor) + &
-                                fse3t(ji,jj,1:k_floor-1)) )
-                  cmax_no=MAXVAL(c_no(2:k_floor))
-
-                  ! number of iterations:
-                  n_iter=min(n_itermax,int(cmax_no)+1)
-                  IF (ln_ctl.AND.(n_iter .gt. 1)) THEN
-                      WRITE(numout,*) 'vertical_movement_fabm():'
-                      WRITE(numout,*) '   Maximum Courant number is ',cmax_no,'.'
-                      WRITE(numout,*) '   ',n_iter,' iterations used for vertical advection.'
-                  ENDIF
-
-                  ! effective Courant number:
-                  c_no=c_no/n_iter
-
-                  tr_it(1:k_floor)=trb(ji,jj,1:k_floor,jp_fabm_m1+jn)
-                  DO n_count=1,n_iter ! Iterative loop
-                     !Compute slope ratio
-                     IF (k_floor.gt.2) THEN !More than 2 vertical wet points
-                        IF (k_floor.gt.3) THEN
-                          WHERE (w_if(3:k_floor-1,jn).ge.0._wp) !upward movement
-                           tr_u(3:k_floor-1)=tr_it(4:k_floor)
-                           tr_c(3:k_floor-1)=tr_it(3:k_floor-1)
-                           tr_d(3:k_floor-1)=tr_it(2:k_floor-2)
-                          ELSEWHERE !downward movement
-                           tr_u(3:k_floor-1)=tr_it(1:k_floor-3)
-                           tr_c(3:k_floor-1)=tr_it(2:k_floor-2)
-                           tr_d(3:k_floor-1)=tr_it(3:k_floor-1)
-                          ENDWHERE
-                        ENDIF
-                        IF (w_if(2,jn).ge.0._wp) THEN
-                           tr_u(2)=tr_it(3)
-                           tr_c(2)=tr_it(2)
-                           tr_d(2)=tr_it(1)
-                        ELSE
-                           tr_u(2)=tr_it(1)
-                           tr_c(2)=tr_it(1)
-                           tr_d(2)=tr_it(2)
-                        ENDIF
-                        IF (w_if(k_floor,jn).ge.0._wp) THEN
-                           tr_u(k_floor)=tr_it(k_floor)
-                           tr_c(k_floor)=tr_it(k_floor)
-                           tr_d(k_floor)=tr_it(k_floor-1)
-                        ELSE
-                           tr_u(k_floor)=tr_it(k_floor-2)
-                           tr_c(k_floor)=tr_it(k_floor-1)
-                           tr_d(k_floor)=tr_it(k_floor)
-                        ENDIF
-                     ELSE !only 2 vertical wet points, i.e. only 1 interface
-                        IF (w_if(k_floor,jn).ge.0._wp) THEN
-                           tr_u(2)=tr_it(2)
-                           tr_c(2)=tr_it(2)
-                           tr_d(2)=tr_it(1)
-                        ELSE
-                           tr_u(2)=tr_it(1)
-                           tr_c(2)=tr_it(1)
-                           tr_d(2)=tr_it(2)
-                        ENDIF
-                     ENDIF
-                     WHERE (abs(tr_d(2:k_floor)-tr_c(2:k_floor)).gt.1.e-10_wp)
-                        tr_slope(2:k_floor)= &
-                           (tr_c(2:k_floor)-tr_u(2:k_floor))/ &
-                           (tr_d(2:k_floor)-tr_c(2:k_floor))
-                     ELSEWHERE
-                        tr_slope(2:k_floor)=SIGN(1._wp,w_if(2:k_floor,jn))* &
-                              (tr_c(2:k_floor)-tr_u(2:k_floor))*1.e10_wp
-                     ENDWHERE
-
-                     !QUICKEST flux limiter:
-                     x_fac(2:k_floor)=(1._wp-2._wp*c_no(2:k_floor))/6._wp
-                     phi_lim(2:k_floor)=(0.5_wp+x_fac(2:k_floor)) + &
-                        (0.5_wp-x_Fac(2:k_floor))*tr_slope(2:k_floor)
-                     flux_lim(2:k_floor)=max( 0._wp, &
-                       min( phi_lim(2:k_floor),2._wp/(1._wp-c_no(2:k_floor)), &
-                         2._wp*tr_slope(2:k_floor)/(c_no(2:k_floor)+1.e-10_wp)) )
-
-                     ! Compute limited flux:
-                     flux_if(2:k_floor,jn) = w_if(2:k_floor,jn)* &
-                        ( tr_c(2:k_floor) + &
-                        0.5_wp*flux_lim(2:k_floor)*(1._wp-c_no(2:k_floor))* &
-                        (tr_d(2:k_floor)-tr_c(2:k_floor)) )
-
-                     ! Compute pseudo update for trend aggregation:
-                     tr_it(1:k_floor-1) = tr_it(1:k_floor-1) + &
-                        z2dt/float(n_iter)/fse3t(ji,jj,1:k_floor-1)* &
-                        flux_if(2:k_floor,jn)
-                     tr_it(2:k_floor) = tr_it(2:k_floor) - &
-                        z2dt/float(n_iter)/fse3t(ji,jj,2:k_floor)* &
-                        flux_if(2:k_floor,jn)
-
-                  ENDDO ! Iterative loop
-
-                  ! Estimate rate of change from pseudo state updates (source
-                  ! splitting):
-                  tra(ji,jj,1:k_floor,jp_fabm_m1+jn) = &
-                     tra(ji,jj,1:k_floor,jp_fabm_m1+jn) + &
-                     (tr_it(1:k_floor) - trb(ji,jj,1:k_floor,jp_fabm_m1+jn))/z2dt
-#if defined key_trdtrc && defined key_iomput
-                  IF( lk_trdtrc .AND. ln_trdtrc( jp_fabm_m1+jn ) ) THEN
-                    tr_vmv(ji,jj,1:k_floor,jn)=(tr_it(1:k_floor) - trb(ji,jj,1:k_floor,jn))/z2dt
+                  IF (ALL(w_ct(ji,1:k_floor,jn) == 0._wp)) CYCLE
+
+                  ! Compute velocities at interfaces
+                  w_if(1:k_floor-1) = zwgt_if(1:k_floor-1) * w_ct(ji,1:k_floor-1,jn) + (1._wp - zwgt_if(1:k_floor-1)) * w_ct(ji,2:k_floor,jn)
+
+                  ! Compute change (per volume) due to vertical movement per layer
+                  IF (method == 1) THEN
+                     CALL advect_1(k_floor, trn(ji,jj,1:k_floor,jp_fabm_m1+jn), w_if(1:k_floor-1), h(1:k_floor), z2dt, dc(1:k_floor))
+                  ELSE
+                     CALL advect_3(k_floor, trb(ji,jj,1:k_floor,jp_fabm_m1+jn), w_if(1:k_floor-1), h(1:k_floor), z2dt, dc(1:k_floor))
                   END IF
-#endif
-               ENDDO ! State loop
+
+                  ! Incorporate change due to vertical movement in sources-sinks
+                  tra(ji,jj,1:k_floor,jp_fabm_m1+jn) = tra(ji,jj,1:k_floor,jp_fabm_m1+jn) + dc(1:k_floor)
+
+#if defined key_trdtrc && defined key_iomput
+                  ! Store change due to vertical movement as diagnostic
+                  IF( lk_trdtrc .AND. ln_trdtrc( jp_fabm_m1+jn)) tr_vmv(ji,jj,1:k_floor,jn) = dc(1:k_floor)
+#endif
+               END DO ! State loop
             END IF ! Level check
          END DO ! i-loop
       END DO ! j-loop
+
 #if defined key_trdtrc && defined key_iomput
       DO jn=1,jp_fabm ! State loop
@@ -220,4 +123,158 @@
    END SUBROUTINE compute_vertical_movement
 
+   SUBROUTINE advect_1(nk, c, w, h, dt, trend)
+      INTEGER,  INTENT(IN)  :: nk
+      REAL(wp), INTENT(IN)  :: c(1:nk)
+      REAL(wp), INTENT(IN)  :: w(1:nk-1)
+      REAL(wp), INTENT(IN)  :: h(1:nk)
+      REAL(wp), INTENT(IN)  :: dt
+      REAL(wp), INTENT(OUT) :: trend(1:nk)
+
+      REAL(wp) :: flux(0:nk)
+      INTEGER  :: jk
+      ! Compute fluxes (per surface area) over at interfaces (remember: positive for upwards)
+      flux(0) = 0._wp
+      DO jk=1,nk-1 ! k-loop
+         IF (w(jk) > 0) THEN
+            ! Upward movement (source layer is jk+1)
+            flux(jk) = min(w(jk), h(jk+1)/dt) * c(jk+1)
+         ELSE
+            ! Downward movement (source layer is jk)
+            flux(jk) = max(w(jk), -h(jk)/dt) * c(jk)
+         END IF
+      END DO
+      flux(nk) = 0._wp
+      trend = (flux(1:nk) - flux(0:nk-1)) / h
+   END SUBROUTINE
+
+   SUBROUTINE advect_3(nk, c_old, w, h, dt, trend)
+      INTEGER,  INTENT(IN)  :: nk
+      REAL(wp), INTENT(IN)  :: c_old(1:nk)
+      REAL(wp), INTENT(IN)  :: w(1:nk-1)
+      REAL(wp), INTENT(IN)  :: h(1:nk)
+      REAL(wp), INTENT(IN)  :: dt
+      REAL(wp), INTENT(OUT) :: trend(1:nk)
+
+      INTEGER, PARAMETER :: n_itermax=100
+      REAL(wp) :: cmax_no
+      REAL(wp) :: cfl(1:nk-1)
+      INTEGER  :: n_iter, n_count, jk
+      REAL(wp) :: c(1:nk)
+      REAL(wp) :: tr_u(1:nk-1)
+      REAL(wp) :: tr_c(1:nk-1)
+      REAL(wp) :: tr_d(1:nk-1)
+      REAL(wp) :: delta_tr_u(1:nk-1)
+      REAL(wp) :: delta_tr(1:nk-1)
+      REAL(wp) :: ratio(1:nk-1)
+      REAL(wp) :: x_fac(1:nk-1)
+      REAL(wp) :: phi_lim(1:nk-1)
+      REAL(wp) :: limiter(1:nk-1)
+      REAL(wp) :: flux_if(1:nk-1)
+
+      c(:) = c_old(:)
+
+      ! get maximum Courant number:
+      cfl = ABS(w) * dt / (0.5_wp * (h(2:nk) + h(1:nk-1)))
+      cmax_no = MAXVAL(cfl)
+
+      ! number of iterations:
+      n_iter = MIN(n_itermax, INT(cmax_no) + 1)
+      IF (ln_ctl.AND.(n_iter .gt. 1)) THEN
+         WRITE(numout,*) 'compute_vertical_movement::advect_3():'
+         WRITE(numout,*) '   Maximum Courant number is ',cmax_no,'.'
+         WRITE(numout,*) '   ',n_iter,' iterations used for vertical advection.'
+      ENDIF
+
+      ! effective Courant number:
+      cfl = cfl/n_iter
+
+      DO n_count=1,n_iter ! Iterative loop
+         ! Determine tracer concentration at 1.5 upstream (tr_u), 0.5 upstream (tr_c), 0.5 downstream (tr_d) from interface
+         IF (nk.gt.2) THEN
+            ! More than 2 vertical wet points
+            IF (nk.gt.3) THEN
+               WHERE (w(2:nk-2).ge.0._wp)
+                  !upward movement
+                  tr_u(2:nk-2)=c(4:nk)
+                  tr_c(2:nk-2)=c(3:nk-1)
+                  tr_d(2:nk-2)=c(2:nk-2)
+               ELSEWHERE
+                  ! downward movement
+                  tr_u(2:nk-2)=c(1:nk-3)
+                  tr_c(2:nk-2)=c(2:nk-2)
+                  tr_d(2:nk-2)=c(3:nk-1)
+               ENDWHERE
+            ENDIF
+
+            ! Interface between surface layer and the next
+            IF (w(1).ge.0._wp) THEN
+               ! upward movement
+               tr_u(1)=c(3)
+               tr_c(1)=c(2)
+               tr_d(1)=c(1)
+            ELSE
+               ! downward movement
+               tr_u(1)=c(1)
+               tr_c(1)=c(1)
+               tr_d(1)=c(2)
+            ENDIF
+
+            ! Interface between bottom layer and the previous
+            IF (w(nk-1).ge.0._wp) THEN
+               ! upward movement
+               tr_u(nk-1)=c(nk)
+               tr_c(nk-1)=c(nk)
+               tr_d(nk-1)=c(nk-1)
+            ELSE
+               ! downward movement
+               tr_u(nk-1)=c(nk-2)
+               tr_c(nk-1)=c(nk-1)
+               tr_d(nk-1)=c(nk)
+            ENDIF
+         ELSE
+            ! only 2 vertical wet points, i.e. only 1 interface
+            IF (w(1).ge.0._wp) THEN
+               ! upward movement
+               tr_u(1)=c(2)
+               tr_c(1)=c(2)
+               tr_d(1)=c(1)
+            ELSE
+               ! downward movement
+               tr_u(1)=c(1)
+               tr_c(1)=c(1)
+               tr_d(1)=c(2)
+            ENDIF
+         ENDIF
+
+         delta_tr_u = tr_c - tr_u
+         delta_tr = tr_d - tr_c
+         WHERE (delta_tr * delta_tr_u > 0._wp)
+            ! Monotonic function over tr_u, tr_c, r_d
+
+            ! Compute slope ratio
+            ratio = delta_tr_u / delta_tr
+
+            ! QUICKEST flux limiter
+            x_fac = (1._wp - 2._wp * cfl) / 6._wp
+            phi_lim = (0.5_wp + x_fac) + (0.5_wp - x_fac) * ratio
+            limiter = MIN(phi_lim, 2._wp / (1._wp - cfl), 2._wp * ratio / (cfl + 1.e-10_wp))
+
+            ! Compute limited flux
+            flux_if = w * (tr_c + 0.5_wp * limiter * (1._wp - cfl) * delta_tr)
+         ELSEWHERE
+            ! Non-monotonic, use 1st order upstream
+            flux_if = w * tr_c
+         ENDWHERE
+
+         ! Compute pseudo update for trend aggregation:
+         c(1:nk-1) = c(1:nk-1) + dt / real(n_iter, wp) / h(1:nk-1) * flux_if
+         c(2:nk)   = c(2:nk)   - dt / real(n_iter, wp) / h(2:nk)   * flux_if
+
+      ENDDO ! Iterative loop
+
+      ! Estimate rate of change from pseudo state updates (source splitting):
+      trend = (c - c_old) / dt
+   END SUBROUTINE
+
 #endif
 END MODULE
