ice_dyn_shared: add optional 'grid_location' argument to seabed_stres…

…s_factor_LKD

In subsequent commits we want to compute the seabed stress factor
'Tb' not only at the U point (NE corner) but also at the E (center of
E face) and N (center of N face) points.

In order to dispatch the computation in this subroutine to different
code paths depending on the grid location (U, N or E), add an optional
argument 'grid_location' that can be used to indicate at which point we
want the factor to be computed. Default it to the 'U' point for
backwards compatibility, such that the existing calls do not have to be
changed.

Note that we need an additional 'l_grid_location' variable to set the
default value, as the Fortran standard does not allow setting the value
of an optional argument if it is not present.

Note also that 'icellu' was incorrectly referencing 'icetmask'. Fix that
by generalizing the description.

cicecore/cicedynB/dynamics/ice_dyn_shared.F90

@@ -930,11 +930,12 @@ subroutine seabed_stress_factor_LKD (nx_block, ny_block,         &
                                            icellu,                     &
                                            indxui,   indxuj,           &
                                            vice,     aice,             &
-                                           hwater,   Tbu)
+                                           hwater,   Tbu,              &
+                                           grid_location)
 
       integer (kind=int_kind), intent(in) :: &
          nx_block, ny_block, &  ! block dimensions
-         icellu                 ! no. of cells where icetmask = 1
+         icellu                 ! no. of cells where ice[uen]mask = 1
 
       integer (kind=int_kind), dimension (nx_block*ny_block), intent(in) :: &
          indxui   , & ! compressed index in i-direction
@@ -946,19 +947,32 @@ subroutine seabed_stress_factor_LKD (nx_block, ny_block,         &
          hwater      ! water depth at tracer location (m)
 
       real (kind=dbl_kind), dimension (nx_block,ny_block), intent(inout) :: &
-         Tbu         ! seabed stress factor (N/m^2)
+         Tbu         ! seabed stress factor at 'grid_location' (N/m^2)
+
+      character(len=*), optional, intent(inout) :: &
+         grid_location    ! grid location (U, E, N), U assumed if not present
 
       real (kind=dbl_kind) :: &
-         au,  & ! concentration of ice at u location
-         hu,  & ! volume per unit area of ice at u location (mean thickness, m)
-         hwu, & ! water depth at u location (m)
-         hcu    ! critical thickness at u location (m)
+         au,  & ! concentration of ice at 'grid_location' 
+         hu,  & ! volume per unit area of ice at 'grid_location' (mean thickness, m)
+         hwu, & ! water depth at 'grid_location' (m)
+         hcu    ! critical thickness at 'grid_location' (m)
 
       integer (kind=int_kind) :: &
          i, j, ij
 
+      character(len=char_len) :: &
+         l_grid_location    ! local version of 'grid_location'
+
       character(len=*), parameter :: subname = '(seabed_stress_factor_LKD)'
 
+      ! Assume U location (NE corner) if grid_location not present
+      if (.not. (present(grid_location))) then
+         l_grid_location = 'U'
+      else
+         l_grid_location = grid_location
+      endif
+
       do ij = 1, icellu
          i = indxui(ij)
          j = indxuj(ij)