diff --git a/libglide/glide_setup.F90 b/libglide/glide_setup.F90 index 6c034925..297dc420 100644 --- a/libglide/glide_setup.F90 +++ b/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 diff --git a/libglide/glide_types.F90 b/libglide/glide_types.F90 index ea81190b..405fa610 100644 --- a/libglide/glide_types.F90 +++ b/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 diff --git a/libglissade/glissade.F90 b/libglissade/glissade.F90 index 562b42d5..64706426 100644 --- a/libglissade/glissade.F90 +++ b/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. diff --git a/libglissade/glissade_bmlt_float.F90 b/libglissade/glissade_bmlt_float.F90 index ba01927f..7479f1ae 100644 --- a/libglissade/glissade_bmlt_float.F90 +++ b/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(&