New option for a uniform thermal forcing anomaly

ISMIP6 ocean forcing is usually applied as a 2D field of thermal forcing.
In experiments to date, thermal forcing anomalies for future warming
have been 2D fields derived from climate models.

This commit adds a simple option to add a uniform thermal forcing anomaly everywhere.
The anomaly is added to the background thermal forcing field.
To apply this option, the user sets thermal_forcing_anomaly to a nonzero value
(in deg C) in the [parameters] section of the config file.

Optionally, the forcing can be phased in linearly by setting the parameter
thermal_forcing_anomaly_timescale to a nonzero value in the config file.
For example, a timescale of 100 years means that after 50 years, only half
the full anomaly is applied.

In the config file, the user should also set the starting time (in years)
for applying the anomaly, if different from 0.  For example, if the anomaly
is applied starting in year 1950, set thermal_forcing_anomaly_tstart = 1950.
(Assuming that the starting time for the anomaly is equal to the time in the
input file does not work for restart runs, since the model tstart variable
is then read from the restart time slice, not the initial input file.)

I tested the new option in 550-year experiments with anomalies of 1 C and 2 C,
to identify the regions most sensitive to a given increase in thermal forcing.

libglide/glide_setup.F90

@@ -2085,6 +2085,9 @@ subroutine handle_parameters(section, model)
 
     ! ocean data parameters
     call GetValue(section, 'gamma0', model%ocean_data%gamma0)
+    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)
 
     ! parameters to adjust input topography
     call GetValue(section, 'adjust_topg_xmin', model%paramets%adjust_topg_xmin)
@@ -2697,6 +2700,14 @@ subroutine print_parameters(model)
     elseif (model%options%whichbmlt_float == BMLT_FLOAT_THERMAL_FORCING) then
        write(message,*) 'gamma0 (nondimensional)  :  ', model%ocean_data%gamma0
        call write_log(message)
+       if (model%ocean_data%thermal_forcing_anomaly /= 0.0d0) then
+          write(message,*) 'thermal forcing anomaly (C) :', model%ocean_data%thermal_forcing_anomaly
+          call write_log(message)
+          write(message,*) 'TF anomaly start time (yr)  :', model%ocean_data%thermal_forcing_anomaly_tstart
+          call write_log(message)
+          write(message,*) 'TF anomaly timescale (yr)   :', model%ocean_data%thermal_forcing_anomaly_timescale
+          call write_log(message)
+       endif
     endif
 
   end subroutine print_parameters

libglide/glide_types.F90

@@ -1742,6 +1742,12 @@ module glide_types
      real(dp), dimension(:,:), pointer :: &
           deltaT_basin => null()                    !> deltaT in each basin (deg C) 
 
+     real(dp) :: &
+          thermal_forcing_anomaly = 0.0d0,  &       !> thermal forcing anomaly (deg C), applied everywhere
+          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.
+
   end type glide_ocean_data
 
   !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

libglissade/glissade.F90

@@ -72,6 +72,7 @@ module glissade
   real(dp), parameter :: thk_init = 500.d0         ! initial thickness (m) for test_transport
   logical, parameter :: test_halo = .false.        ! if true, call test_halo subroutine
 
+  real(dp), parameter :: eps08 = 1.0d-8   ! small number
   real(dp), parameter :: eps11 = 1.0d-11  ! small number
 
 contains
@@ -1420,12 +1421,17 @@ subroutine glissade_bmlt_float_solve(model)
          bmlt_float_transient     ! basal melt rate for ISMIP6 thermal forcing (m/s);
                                   ! take bmlt_float_transient - bmlt_float_baseline to compute anomaly
 
-    real(dp) :: previous_time
+    real(dp) :: previous_time     ! time (yr) at the end of the previous timestep
+    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
+
     integer :: i, j
     integer :: ewn, nsn
     real(dp) :: dew, dns
     integer :: itest, jtest, rtest
 
+
     ! set grid dimensions
     ewn = model%general%ewn
     nsn = model%general%nsn
@@ -1483,6 +1489,30 @@ subroutine glissade_bmlt_float_solve(model)
           print*, 'Compute bmlt_float at runtime from current thermal forcing'
        endif
 
+       !-----------------------------------------------
+       ! Optionally, apply a uniform thermal forcing anomaly everywhere.
+       ! This anomaly can be phased in linearly over a prescribed timescale.
+       !-----------------------------------------------
+
+       if (model%ocean_data%thermal_forcing_anomaly /= 0.0d0) then
+          time_from_start = model%numerics%time - model%ocean_data%thermal_forcing_anomaly_tstart
+          if (time_from_start + eps08 > model%ocean_data%thermal_forcing_anomaly_timescale .or.  &
+               model%ocean_data%thermal_forcing_anomaly_timescale == 0.0d0) then
+             anomaly_fraction = 1.0d0   ! apply the full anomaly
+          else
+             anomaly_fraction = floor(time_from_start + eps08) &
+                  / model%ocean_data%thermal_forcing_anomaly_timescale
+          endif
+          tf_anomaly = anomaly_fraction * model%ocean_data%thermal_forcing_anomaly
+          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
+          endif
+       endif
+
        call glissade_bmlt_float_thermal_forcing(&
             model%options%bmlt_float_thermal_forcing_param, &
             model%options%ocean_data_domain,       &
@@ -1497,7 +1527,8 @@ subroutine glissade_bmlt_float_solve(model)
             model%geometry%lsrf*thk0,              & ! m
             model%geometry%topg*thk0,              & ! m
             model%ocean_data,                      &
-            model%basal_melt%bmlt_float)
+            model%basal_melt%bmlt_float,           &
+            tf_anomaly)                              ! deg C
 
        ! 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

@@ -894,7 +894,8 @@ subroutine glissade_bmlt_float_thermal_forcing(&
        lsrf,                      &
        topg,                      &
        ocean_data,                &
-       bmlt_float)
+       bmlt_float,                &
+       tf_anomaly_in)
 
     use glimmer_paramets, only: thk0, unphys_val
     use glissade_grid_operators, only: glissade_slope_angle
@@ -944,7 +945,10 @@ subroutine glissade_bmlt_float_thermal_forcing(&
                             !> deltaT_basin = temperature corrections per basin for ISMIP6 melt parameterization
 
     real(dp), dimension(:,:), intent(out) :: &
-         bmlt_float                !> basal melt rate for floating ice (m/s)
+         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
 
     ! local variables
 
@@ -958,6 +962,9 @@ subroutine glissade_bmlt_float_thermal_forcing(&
          thermal_forcing_mask,          & ! = 1 where thermal forcing and bmlt_float can be nonzero, else = 0
          new_mask                         ! temporary mask
 
+    real(dp), dimension(ocean_data%nzocn,nx,ny) ::  &
+         thermal_forcing_in               ! TF passed to subroutine interpolate_thermal_forcing_to_lsrf;
+                                          ! optionally corrected for nonzero tf_anomaly
     real(dp), dimension(nx,ny) ::  &
          theta_slope,                   & ! sub-shelf slope angle (radians)
          f_float                          ! weighting function for computing basin averages, in range [0,1]
@@ -967,6 +974,9 @@ subroutine glissade_bmlt_float_thermal_forcing(&
          thermal_forcing_basin,        &  ! basin average thermal forcing (K) at current time
          deltaT_basin_avg                 ! basin average value of deltaT_basin
 
+    real(dp) :: &
+         tf_anomaly                       ! local version of tf_anomaly_in
+
     !TODO - Make H0_float a config parameter?
     real(dp), parameter ::  &
          H0_float = 50.d0                 ! thickness scale (m) for floating ice; used to reduce weights when H < H0_float
@@ -978,6 +988,12 @@ subroutine glissade_bmlt_float_thermal_forcing(&
        print*, '   ocean_data_domain =', ocean_data_domain
     endif
 
+    if (present(tf_anomaly_in)) then
+       tf_anomaly = tf_anomaly_in
+    else
+       tf_anomaly = 0.0d0
+    endif
+
     ! initialize the output
     bmlt_float = 0.0d0
 
@@ -1166,20 +1182,25 @@ subroutine glissade_bmlt_float_thermal_forcing(&
           enddo
        endif
 
-
     endif   ! ocean_data_domain
 
     !-----------------------------------------------
-    ! Interpolate the thermal forcing to the lower ice surface
+    ! Interpolate the thermal forcing to the lower ice surface.
     !-----------------------------------------------
 
+    ! 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
+    endif
+
     call interpolate_thermal_forcing_to_lsrf(&
          nx,                ny,              &
          ocean_data%nzocn,                   &
          ocean_data%zocn,                    &
          thermal_forcing_mask,               &
          lsrf,                               &
-         ocean_data%thermal_forcing,         &
+         thermal_forcing_in,                 &
          ocean_data%thermal_forcing_lsrf)
 
     if (verbose_bmlt_float .and. this_rank==rtest) then