From 2cbb05f90bea396591fd574c4cea51c38a338fb5 Mon Sep 17 00:00:00 2001 From: William Lipscomb Date: Mon, 14 Sep 2020 10:47:51 -0600 Subject: [PATCH] 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 | 11 +++++++++ libglide/glide_types.F90 | 6 +++++ libglissade/glissade.F90 | 35 +++++++++++++++++++++++++++-- libglissade/glissade_bmlt_float.F90 | 31 ++++++++++++++++++++----- 4 files changed, 76 insertions(+), 7 deletions(-) diff --git a/libglide/glide_setup.F90 b/libglide/glide_setup.F90 index 3b9063ad..6c034925 100644 --- a/libglide/glide_setup.F90 +++ b/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 diff --git a/libglide/glide_types.F90 b/libglide/glide_types.F90 index ffa3e9b4..ea81190b 100644 --- a/libglide/glide_types.F90 +++ b/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 !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ diff --git a/libglissade/glissade.F90 b/libglissade/glissade.F90 index c5355cca..562b42d5 100644 --- a/libglissade/glissade.F90 +++ b/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. diff --git a/libglissade/glissade_bmlt_float.F90 b/libglissade/glissade_bmlt_float.F90 index be28b062..ba01927f 100644 --- a/libglissade/glissade_bmlt_float.F90 +++ b/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