Added an option to raise the bed topography in troughs

An earlier commit added the config option 'adjust_input_topography', which could be used
to raise the topography in ocean ridges (i.e., pinning points) while leaving lower
topography unchanged, within a specified rectangular region.

This commit expands that option.  It can now be used to raise the topography in troughs
while leaving higher topography unchanged.

The code decides between these two kinds of adjustment based on the values of the config parameters
topg_max_adjust and topg_no_adjust.  If topg_max_adjust > topg_no_adjust, the high topography is raised.
If topg_max_adjust < topg_no_adjust, the low topography is raised.
(More precisely, the topography is changed by adjust_topg_delta. Positive delta corresponds to raising,
and negative delta to lowering.)

I used this option to raise the topography on the standard 2-km Antarctic grid in the region
with x in the range (513000, 522000) and y in the range (-1659000, -1652000).
In this region there is a low trough which can create a bottleneck and lead to numerical instability.
Raising the trough makes the flow more stable, allowing a longer timestep for ISMIP6-type runs.
For the 2-km grid, the run is stable with 15 steps per year, compared to > 20 steps per year
without the topographic adjustment.

This commit is BFB except when adjust_input_topography = T.

libglide/glide_setup.F90

@@ -2069,8 +2069,8 @@ subroutine handle_parameters(section, model)
     call GetValue(section, 'adjust_topg_xmax', model%paramets%adjust_topg_xmax)
     call GetValue(section, 'adjust_topg_ymin', model%paramets%adjust_topg_ymin)
     call GetValue(section, 'adjust_topg_ymax', model%paramets%adjust_topg_ymax)
-    call GetValue(section, 'adjust_topg_lo',   model%paramets%adjust_topg_lo)
-    call GetValue(section, 'adjust_topg_hi',   model%paramets%adjust_topg_hi)
+    call GetValue(section, 'adjust_topg_no_adjust',  model%paramets%adjust_topg_no_adjust)
+    call GetValue(section, 'adjust_topg_max_adjust', model%paramets%adjust_topg_max_adjust)
     call GetValue(section, 'adjust_topg_delta',   model%paramets%adjust_topg_delta)
 
     ! basal inversion parameters
@@ -2454,9 +2454,9 @@ subroutine print_parameters(model)
        call write_log(message)
        write(message,*) 'adjust_topg_ymax (m)                         : ', model%paramets%adjust_topg_ymax
        call write_log(message)
-       write(message,*) 'adjust_topg_lo (m)                           : ', model%paramets%adjust_topg_lo
+       write(message,*) 'adjust_topg_no_adjust (m)                    : ', model%paramets%adjust_topg_no_adjust
        call write_log(message)
-       write(message,*) 'adjust_topg_hi (m)                           : ', model%paramets%adjust_topg_hi
+       write(message,*) 'adjust_topg_max_adjust (m)                   : ', model%paramets%adjust_topg_max_adjust
        call write_log(message)
        write(message,*) 'adjust_topg_delta (m)                        : ', model%paramets%adjust_topg_delta
        call write_log(message)

libglide/glide_types.F90

@@ -2134,12 +2134,15 @@ module glide_types
                                        !       in regions of rough coastal topography
 
     ! parameters for the adjust_input_topography option
+    ! Note: Depending on the sign of (adjust_topg_max_adjust - adjust_topg_no_adjust), we will either
+    !       (1) change high topography and leave low topography unchanged, or
+    !       (2) change low topography and leave high topography unchanged
     real(dp) :: adjust_topg_xmin = 0.0d0   ! min x value (m) of the region with adjusted topography
     real(dp) :: adjust_topg_xmax = 0.0d0   ! max x value (m) of the region with adjusted topography
     real(dp) :: adjust_topg_ymin = 0.0d0   ! min y value (m) of the region with adjusted topography
     real(dp) :: adjust_topg_ymax = 0.0d0   ! max y value (m) of the region with adjusted topography
-    real(dp) :: adjust_topg_lo = 0.0d0     ! lower threshold (m); adjust where topg > topg_lo
-    real(dp) :: adjust_topg_hi = 0.0d0     ! upper threshold (m); apply max adjustment topg_delta where topg > topg_hi
+    real(dp) :: adjust_topg_no_adjust  = 0.0d0  ! elevation (m) beyond which there is no adjustment of topg
+    real(dp) :: adjust_topg_max_adjust = 0.0d0  ! elevation (m) beyond which there is full adjustment of topg
     real(dp) :: adjust_topg_delta = 0.0d0  ! max adjustment, applied where topg > topg_hi
 
   end type glide_paramets

libglissade/glissade_utils.F90

@@ -77,7 +77,7 @@ subroutine glissade_adjust_thickness(model)
     integer :: nx, ny
     integer :: itest, jtest, rtest
 
-    logical, parameter :: verbose_adjust_thickness = .false.
+    logical, parameter :: verbose_adjust_thickness = .true.
 
     ! Copy some model variables to local variables
 
@@ -390,8 +390,21 @@ subroutine glissade_adjust_topography(model)
     integer :: itest, jtest, rtest
 
     real(dp) :: factor
-    real(dp) :: xmin, xmax, ymin, ymax
-    real(dp) :: topg_lo, topg_hi, topg_delta
+
+    real(dp) :: &
+         xmin, xmax, ymin, ymax  ! x and y boundaries of adjusted region
+
+    ! Note: There are two ways to use these parameters.
+    ! (1) If topg_max_adjust > topg_no_adjust, then we change high topography (topg > topg_max_adjust)
+    !     by topg_delta and leave low topography (topg < topg_no_adjust) unchanged.
+    ! (2) If topg_max_adjust < topg_no_adjust, then we change low topography (topg < topg_max_adjust)
+    !     by topg_delta and leave high topography (topg > topg_no_adjust) unchanged.
+    ! Between topg_no_adjust and topg_max_adjust, the adjustment is phased in linearly.
+
+    real(dp) :: &
+         topg_no_adjust, &    ! elevation (m) beyond which there is no adjustment
+         topg_max_adjust, &   ! elevation (m) beyond which there is full adjustment (by topg_delta)
+         topg_delta           ! max change in topography (m); can be either sign
 
     logical, parameter :: verbose_adjust_topg = .true.
 
@@ -413,8 +426,8 @@ subroutine glissade_adjust_topography(model)
     xmax = model%paramets%adjust_topg_xmax
     ymin = model%paramets%adjust_topg_ymin
     ymax = model%paramets%adjust_topg_ymax
-    topg_lo = model%paramets%adjust_topg_lo
-    topg_hi = model%paramets%adjust_topg_hi
+    topg_no_adjust = model%paramets%adjust_topg_no_adjust
+    topg_max_adjust = model%paramets%adjust_topg_max_adjust
     topg_delta = model%paramets%adjust_topg_delta
 
     if (verbose_adjust_topg .and. this_rank == rtest) then
@@ -426,7 +439,7 @@ subroutine glissade_adjust_topography(model)
        print*, 'thck, topg =', model%geometry%thck(i,j)*thk0, model%geometry%topg(i,j)*thk0
        print*, 'xmin, xmax =', xmin, xmax
        print*, 'ymin, ymax =', ymin, ymax
-       print*, 'topg_lo, topg_hi (m) =', topg_lo, topg_hi
+       print*, 'topg_no_adjust, topg_max_adjust (m) =', topg_no_adjust, topg_max_adjust
        print*, 'topg_delta =', topg_delta
     endif
 
@@ -490,21 +503,46 @@ subroutine glissade_adjust_topography(model)
     topg = model%geometry%topg * thk0
 
     ! Apply the topographic correction.
-    ! Where topg > topg_hi, apply the max correction, topg_delta.
-    ! Where topg < topg_lo, apply no correction.
-    ! Where topg_lo < topg < topg_hi, phase in the correction linearly.
-
-    do j = 1, ny
-       do i = 1, nx
-          if (model%general%x1(i) >= xmin .and. model%general%x1(i) <= xmax .and. &
-              model%general%y1(j) >= ymin .and. model%general%y1(j) <= ymax) then
-             if (topg(i,j) > topg_lo) then
-                factor = min((topg(i,j) - topg_lo)/(topg_hi - topg_lo), 1.0d0)
-                topg(i,j) = topg(i,j) + factor * topg_delta
+    ! Case 1: topg_max_adjust > topg_no_adjust; change higher topography and leave lower topography unchanged
+    ! Case 2: topg_max_adjust < topg_no_adjust; change lower topography and leave higher topography unchanged
+
+    if (topg_max_adjust > topg_no_adjust) then
+
+       ! Where topg > topg_max_adjust, apply the max correction, topg_delta.
+       ! Where topg < topg_no_adjust, apply no correction.
+       ! Where topg_no_adjust < topg < topg_max_adjust, phase in the correction linearly.
+
+       do j = 1, ny
+          do i = 1, nx
+             if (model%general%x1(i) >= xmin .and. model%general%x1(i) <= xmax .and. &
+                 model%general%y1(j) >= ymin .and. model%general%y1(j) <= ymax) then
+                if (topg(i,j) > topg_no_adjust) then
+                   factor = min((topg(i,j) - topg_no_adjust)/(topg_max_adjust - topg_no_adjust), 1.0d0)
+                   topg(i,j) = topg(i,j) + factor * topg_delta
+                endif
+             endif
+          enddo
+       enddo
+
+    elseif (topg_max_adjust < topg_no_adjust) then
+
+       ! Where topg < topg_max_adjust, apply the max correction, topg_delta.
+       ! Where topg > topg_no_adjust, apply no correction.
+       ! Where topg_max_adjust < topg < topg_no_adjust, phase in the correction linearly.
+
+       do j = 1, ny
+          do i = 1, nx
+             if (model%general%x1(i) >= xmin .and. model%general%x1(i) <= xmax .and. &
+                 model%general%y1(j) >= ymin .and. model%general%y1(j) <= ymax) then
+                if (topg(i,j) < topg_no_adjust) then
+                   factor = min((topg_no_adjust - topg(i,j))/(topg_no_adjust - topg_max_adjust), 1.0d0)
+                   topg(i,j) = topg(i,j) + factor * topg_delta
+                endif
              endif
-          endif
+          enddo
        enddo
-    enddo
+
+    endif
 
     model%geometry%topg = topg / thk0
     deallocate(topg)