From 0f9f48b9eb8a23b29a1adbdf90dba4f5529339de Mon Sep 17 00:00:00 2001
From: Philippe Blain <philippe.blain@canada.ca>
Date: Wed, 17 Nov 2021 12:23:36 -0500
Subject: [PATCH] ice_dyn_shared: add optional 'grid_location' argument to
 seabed_stress_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 | 28 ++++++++++++++-----
 1 file changed, 21 insertions(+), 7 deletions(-)

diff --git a/cicecore/cicedynB/dynamics/ice_dyn_shared.F90 b/cicecore/cicedynB/dynamics/ice_dyn_shared.F90
index 63d116376..85f9a42b9 100755
--- a/cicecore/cicedynB/dynamics/ice_dyn_shared.F90
+++ b/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)