Merge branch 'rescale_pressure' of git://github.com/adcroft/MOM6 into…

… Adcroft_rescale_pressure

  Merge branch 'rescale_pressure' of git://github.com/adcroft/MOM6.  This merge
compiles, but it fails the dimensional consistency testing and will be fixed
shortly.

src/ALE/MOM_regridding.F90

@@ -1386,7 +1386,7 @@ subroutine build_rho_grid( G, GV, US, h, tv, dzInterface, remapCS, CS )
       ! Local depth (G%bathyT is positive)
       nominalDepth = G%bathyT(i,j)*GV%Z_to_H
 
-      call build_rho_column(CS%rho_CS, US, nz, nominalDepth, h(i, j, :), &
+      call build_rho_column(CS%rho_CS, nz, nominalDepth, h(i, j, :), &
                             tv%T(i, j, :), tv%S(i, j, :), tv%eqn_of_state, zNew, &
                             h_neglect=h_neglect, h_neglect_edge=h_neglect_edge)
 
@@ -1501,7 +1501,7 @@ subroutine build_grid_HyCOM1( G, GV, US, h, tv, h_new, dzInterface, CS )
              ( 0.5 * ( z_col(K) + z_col(K+1) ) * (GV%H_to_RZ*GV%g_Earth) - tv%P_Ref )
       enddo
 
-      call build_hycom1_column(CS%hycom_CS, US, tv%eqn_of_state, GV%ke, depth, &
+      call build_hycom1_column(CS%hycom_CS, tv%eqn_of_state, GV%ke, depth, &
                                h(i,j,:), tv%T(i,j,:), tv%S(i,j,:), p_col, &
                                z_col, z_col_new, zScale=GV%Z_to_H, &
                                h_neglect=h_neglect, h_neglect_edge=h_neglect_edge)
@@ -1636,7 +1636,7 @@ subroutine build_grid_SLight(G, GV, US, h, tv, dzInterface, CS)
                     ( 0.5 * ( z_col(K) + z_col(K+1) ) * (GV%H_to_RZ*GV%g_Earth) - tv%P_Ref )
       enddo
 
-      call build_slight_column(CS%slight_CS, US, tv%eqn_of_state, GV%H_to_RZ*GV%g_Earth, &
+      call build_slight_column(CS%slight_CS, tv%eqn_of_state, GV%H_to_RZ*GV%g_Earth, &
                           GV%H_subroundoff, nz, depth, h(i, j, :), &
                           tv%T(i, j, :), tv%S(i, j, :), p_col, z_col, z_col_new, &
                           h_neglect=h_neglect, h_neglect_edge=h_neglect_edge)
@@ -1890,7 +1890,7 @@ subroutine convective_adjustment(G, GV, h, tv)
   do j = G%jsc-1,G%jec+1 ; do i = G%isc-1,G%iec+1
 
     ! Compute densities within current water column
-    call calculate_density( tv%T(i,j,:), tv%S(i,j,:), p_col, densities, tv%eqn_of_state, US=G%US)
+    call calculate_density( tv%T(i,j,:), tv%S(i,j,:), p_col, densities, tv%eqn_of_state)
 
     ! Repeat restratification until complete
     do
@@ -1909,9 +1909,9 @@ subroutine convective_adjustment(G, GV, h, tv)
           tv%S(i,j,k) = S1 ; tv%S(i,j,k+1) = S0
           h(i,j,k)    = h1 ; h(i,j,k+1)    = h0
           ! Recompute densities at levels k and k+1
-          call calculate_density( tv%T(i,j,k), tv%S(i,j,k), p_col(k), densities(k), tv%eqn_of_state, US=G%US)
+          call calculate_density( tv%T(i,j,k), tv%S(i,j,k), p_col(k), densities(k), tv%eqn_of_state)
           call calculate_density( tv%T(i,j,k+1), tv%S(i,j,k+1), p_col(k+1), &
-                                  densities(k+1), tv%eqn_of_state, US=G%US )
+                                  densities(k+1), tv%eqn_of_state )
           stratified = .false.
         endif
       enddo  ! k

src/ALE/coord_adapt.F90

@@ -159,7 +159,7 @@ subroutine build_adapt_column(CS, G, GV, US, tv, i, j, zInt, tInt, sInt, h, zNex
          0.5 * (tInt(i,j,2:nz) + tInt(i,j-1,2:nz)), &
          0.5 * (sInt(i,j,2:nz) + sInt(i,j-1,2:nz)), &
          0.5 * (zInt(i,j,2:nz) + zInt(i,j-1,2:nz)) * (GV%H_to_RZ * GV%g_Earth), &
-         alpha, beta, tv%eqn_of_state, US, dom=(/2,nz/))
+         alpha, beta, tv%eqn_of_state, (/2,nz/))
 
     del2sigma(2:nz) = del2sigma(2:nz) + &
          (alpha(2:nz) * (tInt(i,j-1,2:nz) - tInt(i,j,2:nz)) + &
@@ -171,7 +171,7 @@ subroutine build_adapt_column(CS, G, GV, US, tv, i, j, zInt, tInt, sInt, h, zNex
          0.5 * (tInt(i,j,2:nz) + tInt(i,j+1,2:nz)), &
          0.5 * (sInt(i,j,2:nz) + sInt(i,j+1,2:nz)), &
          0.5 * (zInt(i,j,2:nz) + zInt(i,j+1,2:nz)) * (GV%H_to_RZ * GV%g_Earth), &
-         alpha, beta, tv%eqn_of_state, US, dom=(/2,nz/))
+         alpha, beta, tv%eqn_of_state, (/2,nz/))
 
     del2sigma(2:nz) = del2sigma(2:nz) + &
          (alpha(2:nz) * (tInt(i,j+1,2:nz) - tInt(i,j,2:nz)) + &
@@ -183,7 +183,7 @@ subroutine build_adapt_column(CS, G, GV, US, tv, i, j, zInt, tInt, sInt, h, zNex
          0.5 * (tInt(i,j,2:nz) + tInt(i-1,j,2:nz)), &
          0.5 * (sInt(i,j,2:nz) + sInt(i-1,j,2:nz)), &
          0.5 * (zInt(i,j,2:nz) + zInt(i-1,j,2:nz)) * (GV%H_to_RZ * GV%g_Earth), &
-         alpha, beta, tv%eqn_of_state, US, dom=(/2,nz/))
+         alpha, beta, tv%eqn_of_state, (/2,nz/))
 
     del2sigma(2:nz) = del2sigma(2:nz) + &
          (alpha(2:nz) * (tInt(i-1,j,2:nz) - tInt(i,j,2:nz)) + &
@@ -195,7 +195,7 @@ subroutine build_adapt_column(CS, G, GV, US, tv, i, j, zInt, tInt, sInt, h, zNex
          0.5 * (tInt(i,j,2:nz) + tInt(i+1,j,2:nz)), &
          0.5 * (sInt(i,j,2:nz) + sInt(i+1,j,2:nz)), &
          0.5 * (zInt(i,j,2:nz) + zInt(i+1,j,2:nz)) * (GV%H_to_RZ * GV%g_Earth), &
-         alpha, beta, tv%eqn_of_state, US, dom=(/2,nz/))
+         alpha, beta, tv%eqn_of_state, (/2,nz/))
 
     del2sigma(2:nz) = del2sigma(2:nz) + &
          (alpha(2:nz) * (tInt(i+1,j,2:nz) - tInt(i,j,2:nz)) + &
@@ -209,7 +209,7 @@ subroutine build_adapt_column(CS, G, GV, US, tv, i, j, zInt, tInt, sInt, h, zNex
   ! a positive curvature means we're too light relative to adjacent columns,
   ! so del2sigma needs to be positive too (push the interface deeper)
   call calculate_density_derivs(tInt(i,j,:), sInt(i,j,:), zInt(i,j,:) * (GV%H_to_RZ * GV%g_Earth), &
-       alpha, beta, tv%eqn_of_state, US, dom=(/1,nz/)) !### This should be (/1,nz+1/) - see 25 lines below.
+       alpha, beta, tv%eqn_of_state, (/1,nz/)) !### This should be (/1,nz+1/) - see 25 lines below.
   do K = 2, nz
     ! TODO make lower bound here configurable
     dh_d2s(K) = del2sigma(K) * (0.5 * (h(i,j,k-1) + h(i,j,k))) / &

src/ALE/coord_hycom.F90

@@ -4,7 +4,6 @@ module coord_hycom
 ! This file is part of MOM6. See LICENSE.md for the license.
 
 use MOM_error_handler, only : MOM_error, FATAL
-use MOM_unit_scaling,  only : unit_scale_type
 use MOM_EOS,           only : EOS_type, calculate_density
 use regrid_interp,     only : interp_CS_type, build_and_interpolate_grid
 
@@ -96,10 +95,9 @@ subroutine set_hycom_params(CS, max_interface_depths, max_layer_thickness, inter
 end subroutine set_hycom_params
 
 !> Build a HyCOM coordinate column
-subroutine build_hycom1_column(CS, US, eqn_of_state, nz, depth, h, T, S, p_col, &
+subroutine build_hycom1_column(CS, eqn_of_state, nz, depth, h, T, S, p_col, &
                                z_col, z_col_new, zScale, h_neglect, h_neglect_edge)
   type(hycom_CS),        intent(in)    :: CS    !< Coordinate control structure
-  type(unit_scale_type), intent(in)    :: US    !< A dimensional unit scaling type
   type(EOS_type),        pointer       :: eqn_of_state !< Equation of state structure
   integer,               intent(in)    :: nz    !< Number of levels
   real,                  intent(in)    :: depth !< Depth of ocean bottom (positive [H ~> m or kg m-2])
@@ -133,7 +131,7 @@ subroutine build_hycom1_column(CS, US, eqn_of_state, nz, depth, h, T, S, p_col,
   z_scale = 1.0 ; if (present(zScale)) z_scale = zScale
 
   ! Work bottom recording potential density
-  call calculate_density(T, S, p_col, rho_col, eqn_of_state, US=US)
+  call calculate_density(T, S, p_col, rho_col, eqn_of_state)
   ! This ensures the potential density profile is monotonic
   ! although not necessarily single valued.
   do k = nz-1, 1, -1

src/ALE/coord_rho.F90

@@ -5,7 +5,6 @@ module coord_rho
 
 use MOM_error_handler, only : MOM_error, FATAL
 use MOM_remapping,     only : remapping_CS, remapping_core_h
-use MOM_unit_scaling,  only : unit_scale_type
 use MOM_EOS,           only : EOS_type, calculate_density
 use regrid_interp,     only : interp_CS_type, build_and_interpolate_grid, DEGREE_MAX
 
@@ -88,10 +87,9 @@ end subroutine set_rho_params
 !!
 !! 1. Density profiles are calculated on the source grid.
 !! 2. Positions of target densities (for interfaces) are found by interpolation.
-subroutine build_rho_column(CS, US, nz, depth, h, T, S, eqn_of_state, z_interface, &
+subroutine build_rho_column(CS, nz, depth, h, T, S, eqn_of_state, z_interface, &
                             h_neglect, h_neglect_edge)
   type(rho_CS),        intent(in)    :: CS !< coord_rho control structure
-  type(unit_scale_type), intent(in)  :: US    !< A dimensional unit scaling type
   integer,             intent(in)    :: nz !< Number of levels on source grid (i.e. length of  h, T, S)
   real,                intent(in)    :: depth !< Depth of ocean bottom (positive in m)
   real, dimension(nz), intent(in)    :: h  !< Layer thicknesses [H ~> m or kg m-2]
@@ -126,7 +124,7 @@ subroutine build_rho_column(CS, US, nz, depth, h, T, S, eqn_of_state, z_interfac
 
     ! Compute densities on source column
     pres(:) = CS%ref_pressure
-    call calculate_density(T, S, pres, densities, eqn_of_state, US=US)
+    call calculate_density(T, S, pres, densities, eqn_of_state)
     do k = 1,count_nonzero_layers
       densities(k) = densities(mapping(k))
     enddo
@@ -185,10 +183,9 @@ end subroutine build_rho_column
 !! 4. T & S are remapped onto the new grid.
 !! 5. Return to step 1 until convergence or until the maximum number of
 !!    iterations is reached, whichever comes first.
-subroutine build_rho_column_iteratively(CS, US, remapCS, nz, depth, h, T, S, eqn_of_state, &
+subroutine build_rho_column_iteratively(CS, remapCS, nz, depth, h, T, S, eqn_of_state, &
                                         zInterface, h_neglect, h_neglect_edge, dev_tol)
   type(rho_CS),          intent(in)    :: CS !< Regridding control structure
-  type(unit_scale_type), intent(in)    :: US !< A dimensional unit scaling type
   type(remapping_CS),    intent(in)    :: remapCS !< Remapping parameters and options
   integer,               intent(in)    :: nz !< Number of levels
   real,                  intent(in)    :: depth !< Depth of ocean bottom [Z ~> m]
@@ -250,7 +247,7 @@ subroutine build_rho_column_iteratively(CS, US, remapCS, nz, depth, h, T, S, eqn
     enddo
 
     ! Compute densities within current water column
-    call calculate_density( T_tmp, S_tmp, pres, densities, eqn_of_state, US=US)
+    call calculate_density( T_tmp, S_tmp, pres, densities, eqn_of_state)
 
     do k = 1,count_nonzero_layers
       densities(k) = densities(mapping(k))

src/ALE/coord_slight.F90

@@ -4,7 +4,6 @@ module coord_slight
 ! This file is part of MOM6. See LICENSE.md for the license.
 
 use MOM_error_handler, only : MOM_error, FATAL
-use MOM_unit_scaling,  only : unit_scale_type
 use MOM_EOS,           only : EOS_type, calculate_compress
 use MOM_EOS,           only : calculate_density, calculate_density_derivs
 use regrid_interp,     only : interp_CS_type, regridding_set_ppolys
@@ -178,11 +177,10 @@ subroutine set_slight_params(CS, max_interface_depths, max_layer_thickness, &
 end subroutine set_slight_params
 
 !> Build a SLight coordinate column
-subroutine build_slight_column(CS, US, eqn_of_state, H_to_pres, H_subroundoff, &
+subroutine build_slight_column(CS, eqn_of_state, H_to_pres, H_subroundoff, &
                                nz, depth, h_col, T_col, S_col, p_col, z_col, z_col_new, &
                                h_neglect, h_neglect_edge)
   type(slight_CS),       intent(in)    :: CS    !< Coordinate control structure
-  type(unit_scale_type), intent(in)    :: US    !< A dimensional unit scaling type
   type(EOS_type),        pointer       :: eqn_of_state !< Equation of state structure
   real,                  intent(in)    :: H_to_pres !< A conversion factor from thicknesses to
                                                 !! scaled pressure [R L2 T-2 H-1 ~> Pa m-1 or Pa m2 kg-1]
@@ -253,7 +251,7 @@ subroutine build_slight_column(CS, US, eqn_of_state, H_to_pres, H_subroundoff, &
     dz = (z_col(nz+1) - z_col(1)) / real(nz)
     do K=2,nz ; z_col_new(K) = z_col(1) + dz*real(K-1) ; enddo
   else
-    call calculate_density(T_col, S_col, p_col, rho_col, eqn_of_state, US=US)
+    call calculate_density(T_col, S_col, p_col, rho_col, eqn_of_state)
 
     ! Find the locations of the target potential densities, flagging
     ! locations in apparently unstable regions as not reliable.
@@ -375,11 +373,11 @@ subroutine build_slight_column(CS, US, eqn_of_state, H_to_pres, H_subroundoff, &
       T_int(nz+1) = T_f(nz) ; S_int(nz+1) = S_f(nz)
       p_IS(nz+1) = z_col(nz+1) * H_to_pres
       call calculate_density_derivs(T_int, S_int, p_IS, drhoIS_dT, drhoIS_dS, &
-                                    eqn_of_state, US, dom=(/2,nz/))
+                                    eqn_of_state, dom=(/2,nz/))
       call calculate_density_derivs(T_int, S_int, p_R, drhoR_dT, drhoR_dS, &
-                                    eqn_of_state, US, dom=(/2,nz/))
+                                    eqn_of_state, dom=(/2,nz/))
       if (CS%compressibility_fraction > 0.0) then
-        call calculate_compress(T_int, S_int, p_R(:), rho_tmp, drho_dp, 2, nz-1, eqn_of_state, US)
+        call calculate_compress(T_int, S_int, p_R(:), rho_tmp, drho_dp, 2, nz-1, eqn_of_state)
       else
         do K=2,nz ; drho_dp(K) = 0.0 ; enddo
       endif

src/core/MOM.F90

@@ -2242,7 +2242,7 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, restart_CSp, &
   ! Use the Wright equation of state by default, unless otherwise specified
   ! Note: this line and the following block ought to be in a separate
   ! initialization routine for tv.
-  if (use_EOS) call EOS_init(param_file, CS%tv%eqn_of_state)
+  if (use_EOS) call EOS_init(param_file, CS%tv%eqn_of_state, US)
   if (use_temperature) then
     allocate(CS%tv%TempxPmE(isd:ied,jsd:jed)) ; CS%tv%TempxPmE(:,:) = 0.0
     if (use_geothermal) then
@@ -2905,7 +2905,7 @@ subroutine adjust_ssh_for_p_atm(tv, G, GV, US, ssh, p_atm, use_EOS)
     do j=js,je
       if (calc_rho) then
         call calculate_density(tv%T(:,j,1), tv%S(:,j,1), 0.5*p_atm(:,j), Rho_conv, &
-                               tv%eqn_of_state, US, dom=EOS_domain(G%HI))
+                               tv%eqn_of_state, dom=EOS_domain(G%HI))
         do i=is,ie
           IgR0 = US%Z_to_m / (Rho_conv(i) * GV%g_Earth)
           ssh(i,j) = ssh(i,j) + p_atm(i,j) * IgR0

src/core/MOM_PressureForce_Montgomery.F90

@@ -188,7 +188,7 @@ subroutine PressureForce_Mont_nonBouss(h, tv, PFu, PFv, G, GV, US, CS, p_atm, pb
       !$OMP parallel do default(shared)
       do k=1,nz
         call int_specific_vol_dp(tv%T(:,:,k), tv%S(:,:,k), p(:,:,k), p(:,:,k+1), &
-                                 0.0, G%HI, tv%eqn_of_state, dz_geo(:,:,k), halo_size=1, US=US)
+                                 0.0, G%HI, tv%eqn_of_state, dz_geo(:,:,k), halo_size=1)
       enddo
       !$OMP parallel do default(shared)
       do j=Jsq,Jeq+1 ; do k=1,nz ; do i=Isq,Ieq+1
@@ -230,7 +230,7 @@ subroutine PressureForce_Mont_nonBouss(h, tv, PFu, PFv, G, GV, US, CS, p_atm, pb
           tv_tmp%T(i,j,k) = tv%T(i,j,k) ; tv_tmp%S(i,j,k) = tv%S(i,j,k)
         enddo ; enddo
         call calculate_density(tv%T(:,j,nkmb), tv%S(:,j,nkmb), p_ref, Rho_cv_BL(:), &
-                               tv%eqn_of_state, US=US, dom=EOSdom)
+                               tv%eqn_of_state, dom=EOSdom)
         do k=nkmb+1,nz ; do i=Isq,Ieq+1
           if (GV%Rlay(k) < Rho_cv_BL(i)) then
             tv_tmp%T(i,j,k) = tv%T(i,j,nkmb) ; tv_tmp%S(i,j,k) = tv%S(i,j,nkmb)
@@ -247,7 +247,7 @@ subroutine PressureForce_Mont_nonBouss(h, tv, PFu, PFv, G, GV, US, CS, p_atm, pb
     !$OMP parallel do default(shared) private(rho_in_situ)
     do k=1,nz ; do j=Jsq,Jeq+1
       call calculate_density(tv_tmp%T(:,j,k), tv_tmp%S(:,j,k), p_ref, rho_in_situ, &
-                             tv%eqn_of_state, US=US, dom=EOSdom)
+                             tv%eqn_of_state, dom=EOSdom)
       do i=Isq,Ieq+1 ; alpha_star(i,j,k) = 1.0 / rho_in_situ(i) ; enddo
     enddo ; enddo
   endif                                               ! use_EOS
@@ -487,7 +487,7 @@ subroutine PressureForce_Mont_Bouss(h, tv, PFu, PFv, G, GV, US, CS, p_atm, pbce,
           tv_tmp%T(i,j,k) = tv%T(i,j,k) ; tv_tmp%S(i,j,k) = tv%S(i,j,k)
         enddo ; enddo
         call calculate_density(tv%T(:,j,nkmb), tv%S(:,j,nkmb), p_ref, Rho_cv_BL(:), &
-                               tv%eqn_of_state, US=US, dom=EOSdom)
+                               tv%eqn_of_state, dom=EOSdom)
 
         do k=nkmb+1,nz ; do i=Isq,Ieq+1
           if (GV%Rlay(k) < Rho_cv_BL(i)) then
@@ -508,7 +508,7 @@ subroutine PressureForce_Mont_Bouss(h, tv, PFu, PFv, G, GV, US, CS, p_atm, pbce,
     !$OMP parallel do default(shared)
     do k=1,nz+1 ; do j=Jsq,Jeq+1
       call calculate_density(tv_tmp%T(:,j,k), tv_tmp%S(:,j,k), p_ref, rho_star(:,j,k), &
-                             tv%eqn_of_state, US=US, dom=EOSdom)
+                             tv%eqn_of_state, dom=EOSdom)
       do i=Isq,Ieq+1 ; rho_star(i,j,k) = G_Rho0*rho_star(i,j,k) ; enddo
     enddo ; enddo
   endif                                               ! use_EOS
@@ -666,7 +666,7 @@ subroutine Set_pbce_Bouss(e, tv, G, GV, US, Rho0, GFS_scale, pbce, rho_star)
           press(i) = -Rho0xG*e(i,j,1)
         enddo
         call calculate_density(tv%T(:,j,1), tv%S(:,j,1), press, rho_in_situ, &
-                               tv%eqn_of_state, US=US, dom=EOSdom)
+                               tv%eqn_of_state, dom=EOSdom)
         do i=Isq,Ieq+1
           pbce(i,j,1) = G_Rho0*(GFS_scale * rho_in_situ(i)) * GV%H_to_Z
         enddo
@@ -677,7 +677,7 @@ subroutine Set_pbce_Bouss(e, tv, G, GV, US, Rho0, GFS_scale, pbce, rho_star)
             S_int(i) = 0.5*(tv%S(i,j,k-1)+tv%S(i,j,k))
           enddo
           call calculate_density_derivs(T_int, S_int, press, dR_dT, dR_dS, &
-                                        tv%eqn_of_state, US=US, dom=EOSdom)
+                                        tv%eqn_of_state, dom=EOSdom)
           do i=Isq,Ieq+1
             pbce(i,j,k) = pbce(i,j,k-1) + G_Rho0 * &
                ((e(i,j,K) - e(i,j,nz+1)) * Ihtot(i)) * &
@@ -764,7 +764,7 @@ subroutine Set_pbce_nonBouss(p, tv, G, GV, US, GFS_scale, pbce, alpha_star)
       !$OMP parallel do default(shared) private(T_int,S_int,dR_dT,dR_dS,rho_in_situ)
       do j=Jsq,Jeq+1
         call calculate_density(tv%T(:,j,nz), tv%S(:,j,nz), p(:,j,nz+1), rho_in_situ, &
-                               tv%eqn_of_state, US=US, dom=EOSdom)
+                               tv%eqn_of_state, dom=EOSdom)
         do i=Isq,Ieq+1
           C_htot(i,j) = dP_dH / ((p(i,j,nz+1)-p(i,j,1)) + dp_neglect)
           pbce(i,j,nz) = dP_dH / (rho_in_situ(i))
@@ -774,9 +774,9 @@ subroutine Set_pbce_nonBouss(p, tv, G, GV, US, GFS_scale, pbce, alpha_star)
             T_int(i) = 0.5*(tv%T(i,j,k)+tv%T(i,j,k+1))
             S_int(i) = 0.5*(tv%S(i,j,k)+tv%S(i,j,k+1))
           enddo
-          call calculate_density(T_int, S_int, p(:,j,k+1), rho_in_situ, tv%eqn_of_state, US=US, dom=EOSdom)
+          call calculate_density(T_int, S_int, p(:,j,k+1), rho_in_situ, tv%eqn_of_state, dom=EOSdom)
           call calculate_density_derivs(T_int, S_int, p(:,j,k+1), dR_dT, dR_dS, &
-                                        tv%eqn_of_state, US=US, dom=EOSdom)
+                                        tv%eqn_of_state, dom=EOSdom)
           do i=Isq,Ieq+1
             pbce(i,j,k) = pbce(i,j,k+1) + ((p(i,j,K+1)-p(i,j,1))*C_htot(i,j)) *  &
                 ((dR_dT(i)*(tv%T(i,j,k+1)-tv%T(i,j,k)) + &