Option to apply thermal forcing anomaly to a single basin

With this commit, the user can specify config variable thermal_forcing_anomaly_basin
to apply the anomaly to a single basin, e.g. the Amundsen Sea (ISMIP6 basin #9).

If this variable is not specified, the default value is 0, which means
the anomaly is applied to all basins.

Applying anomalies to a single basin can be useful on multi-century time scales,
e.g. to determine whether ASE retreat is locally forced or is triggered
by retreat in neighboring basins.

libglide/glide_setup.F90

@@ -2088,6 +2088,7 @@ subroutine handle_parameters(section, model)
     call GetVAlue(section, 'thermal_forcing_anomaly', model%ocean_data%thermal_forcing_anomaly)
     call GetVAlue(section, 'thermal_forcing_anomaly_tstart', model%ocean_data%thermal_forcing_anomaly_tstart)
     call GetVAlue(section, 'thermal_forcing_anomaly_timescale', model%ocean_data%thermal_forcing_anomaly_timescale)
+    call GetVAlue(section, 'thermal_forcing_anomaly_basin', model%ocean_data%thermal_forcing_anomaly_basin)
 
     ! parameters to adjust input topography
     call GetValue(section, 'adjust_topg_xmin', model%paramets%adjust_topg_xmin)
@@ -2707,6 +2708,13 @@ subroutine print_parameters(model)
           call write_log(message)
           write(message,*) 'TF anomaly timescale (yr)   :', model%ocean_data%thermal_forcing_anomaly_timescale
           call write_log(message)
+          if (model%ocean_data%thermal_forcing_anomaly_basin == 0) then
+             call write_log('TF anomaly will be applied to all basins')
+          else
+             write(message,*) 'TF anomaly will be applied to basin', &
+                  model%ocean_data%thermal_forcing_anomaly_basin
+             call write_log(message)
+          endif
        endif
     endif
 

libglide/glide_types.F90

@@ -1747,6 +1747,9 @@ module glide_types
           thermal_forcing_anomaly_tstart = 0.0d0, & !> starting time (yr) for applying or phasing in the anomaly
           thermal_forcing_anomaly_timescale = 0.0d0 !> number of years over which the anomaly is phased in linearly
                                                     !> If set to zero, then the full anomaly is applied immediately.
+     integer :: &
+          thermal_forcing_anomaly_basin = 0         !> basin where anomaly is applied;
+                                                    !> for default value of 0, apply to all basins
 
   end type glide_ocean_data
 

libglissade/glissade.F90

@@ -1425,6 +1425,8 @@ subroutine glissade_bmlt_float_solve(model)
     real(dp) :: time_from_start   ! time (yr) since the start of applying the anomaly
     real(dp) :: anomaly_fraction  ! fraction of full anomaly to apply
     real(dp) :: tf_anomaly        ! uniform thermal forcing anomaly (deg C), applied everywhere
+    real(dp) :: tf_anomaly_basin  ! basin number where anomaly is applied;
+                                  ! for default value of 0, apply to all basins
 
     integer :: i, j
     integer :: ewn, nsn
@@ -1504,12 +1506,16 @@ subroutine glissade_bmlt_float_solve(model)
                   / model%ocean_data%thermal_forcing_anomaly_timescale
           endif
           tf_anomaly = anomaly_fraction * model%ocean_data%thermal_forcing_anomaly
+          tf_anomaly_basin = model%ocean_data%thermal_forcing_anomaly_basin
           if (this_rank == rtest .and. verbose_bmlt_float) then
              print*, 'time_from_start (yr):', time_from_start
              print*, 'thermal forcing anomaly  (deg):', model%ocean_data%thermal_forcing_anomaly
              print*, 'timescale (yr):', model%ocean_data%thermal_forcing_anomaly_timescale
              print*, 'fraction:', anomaly_fraction
              print*, 'current TF anomaly (deg):', tf_anomaly
+             if (model%ocean_data%thermal_forcing_anomaly_timescale /= 0) then
+                print*, 'anomaly applied to basin', model%ocean_data%thermal_forcing_anomaly_basin
+             endif
           endif
        endif
 
@@ -1528,7 +1534,8 @@ subroutine glissade_bmlt_float_solve(model)
             model%geometry%topg*thk0,              & ! m
             model%ocean_data,                      &
             model%basal_melt%bmlt_float,           &
-            tf_anomaly)                              ! deg C
+            tf_anomaly,                            & ! deg C
+            tf_anomaly_basin)
 
        ! There are two ways to compute the transient basal melting from the thermal forcing at runtime:
        ! (1) Use the value just computed, based on the current thermal_forcing.

libglissade/glissade_bmlt_float.F90

@@ -895,7 +895,8 @@ subroutine glissade_bmlt_float_thermal_forcing(&
        topg,                      &
        ocean_data,                &
        bmlt_float,                &
-       tf_anomaly_in)
+       tf_anomaly_in,             &
+       tf_anomaly_basin_in)
 
     use glimmer_paramets, only: thk0, unphys_val
     use glissade_grid_operators, only: glissade_slope_angle
@@ -948,7 +949,8 @@ subroutine glissade_bmlt_float_thermal_forcing(&
          bmlt_float         !> basal melt rate for floating ice (m/s)
 
     real(dp), intent(in), optional ::  &
-         tf_anomaly_in      !> uniform thermal forcing anomaly (deg C), applied everywhere
+         tf_anomaly_in,   &   !> uniform thermal forcing anomaly (deg C), applied everywhere
+         tf_anomaly_basin_in  !> basin where anomaly is applied; for default value of 0, apply to all basins
 
     ! local variables
 
@@ -975,7 +977,8 @@ subroutine glissade_bmlt_float_thermal_forcing(&
          deltaT_basin_avg                 ! basin average value of deltaT_basin
 
     real(dp) :: &
-         tf_anomaly                       ! local version of tf_anomaly_in
+         tf_anomaly,                   &  ! local version of tf_anomaly_in
+         tf_anomaly_basin                 ! local version of tf_anomaly_basin_in
 
     !TODO - Make H0_float a config parameter?
     real(dp), parameter ::  &
@@ -994,6 +997,12 @@ subroutine glissade_bmlt_float_thermal_forcing(&
        tf_anomaly = 0.0d0
     endif
 
+    if (present(tf_anomaly_basin_in)) then
+       tf_anomaly_basin = tf_anomaly_basin_in
+    else
+       tf_anomaly_basin = 0
+    endif
+
     ! initialize the output
     bmlt_float = 0.0d0
 
@@ -1191,7 +1200,17 @@ subroutine glissade_bmlt_float_thermal_forcing(&
     ! Optionally, add a uniform anomaly (= 0 by default) to the thermal forcing.
     thermal_forcing_in = ocean_data%thermal_forcing
     if (tf_anomaly /= 0.0d0) then
-       thermal_forcing_in = ocean_data%thermal_forcing + tf_anomaly
+       if (tf_anomaly_basin >= 1 .and. tf_anomaly_basin <= ocean_data%nbasin) then  ! apply to one basin
+          do j = 1, ny
+             do i = 1, nx
+                if (ocean_data%basin_number(i,j) == tf_anomaly_basin) then
+                   thermal_forcing_in(:,i,j) = thermal_forcing_in(:,i,j) + tf_anomaly
+                endif
+             enddo
+          enddo
+       else   ! apply to all basins
+          thermal_forcing_in = thermal_forcing_in + tf_anomaly
+       endif
     endif
 
     call interpolate_thermal_forcing_to_lsrf(&