From bb8e25b3ec42891ec7385c0cd0897648a3978cfa Mon Sep 17 00:00:00 2001 From: hgoelzer Date: Thu, 6 Feb 2020 11:05:17 +0100 Subject: [PATCH] Adding forced topography options whichrelaxed = 3 and whichrelaxed = 4. They may be used to relax the model to a given topography. --- libglide/glide_setup.F90 | 4 +++ libglide/glide_types.F90 | 10 +++++-- libglide/isostasy.F90 | 60 ++++++++++++++++++++++++++++++++++++++++ libglissade/glissade.F90 | 37 ++++++++++++++++++++++--- 4 files changed, 104 insertions(+), 7 deletions(-) diff --git a/libglide/glide_setup.F90 b/libglide/glide_setup.F90 index 7acbfd24..2e5afb98 100644 --- a/libglide/glide_setup.F90 +++ b/libglide/glide_setup.F90 @@ -2323,6 +2323,10 @@ subroutine print_isostasy(model) call write_log('setting relx to first slice of input topg') elseif (model%isostasy%whichrelaxed==RELAXED_TOPO_COMPUTE) then call write_log('computing relx, given that input topg is in equilibrium') + elseif (model%isostasy%whichrelaxed==RELAXED_TOPO_TARGET) then + call write_log('reading relx as load-independent target topg') + elseif (model%isostasy%whichrelaxed==RELAXED_TOPO_FORCED) then + call write_log('setting relx as load-independent target to first slice of input topg ') else call write_log('Error, unknown whichrelaxed option',GM_FATAL) end if diff --git a/libglide/glide_types.F90 b/libglide/glide_types.F90 index 417e79ce..58905935 100644 --- a/libglide/glide_types.F90 +++ b/libglide/glide_types.F90 @@ -139,10 +139,14 @@ module glide_types integer, parameter :: GTHF_PRESCRIBED_2D = 1 integer, parameter :: GTHF_COMPUTE = 2 - integer, parameter :: RELAXED_TOPO_DEFAULT = 0 ! topo and relx are separate input fields + integer, parameter :: RELAXED_TOPO_DEFAULT = 0 ! topg and relx are separate input fields integer, parameter :: RELAXED_TOPO_INPUT = 1 ! set relx to input topg - integer, parameter :: RELAXED_TOPO_COMPUTE = 2 ! Input topo in isostatic equilibrium - ! compute relaxed topo + integer, parameter :: RELAXED_TOPO_COMPUTE = 2 ! Input topg in isostatic equilibrium + ! compute relx + !HG: Adding options for forced bedrock adjustment + integer, parameter :: RELAXED_TOPO_TARGET = 3 ! use relx as load-independent target for topg + integer, parameter :: RELAXED_TOPO_FORCED = 4 ! set relx to input topg as load-independent + ! target for topg integer, parameter :: ISOSTASY_NONE = 0 integer, parameter :: ISOSTASY_COMPUTE = 1 diff --git a/libglide/isostasy.F90 b/libglide/isostasy.F90 index e7c811bd..67d2e152 100644 --- a/libglide/isostasy.F90 +++ b/libglide/isostasy.F90 @@ -86,6 +86,24 @@ module isostasy ! ! In general, the preferred setting is whichrelaxed = 0, with topg and relx read in separately ! from the input file. The other settings have specialized uses but may be inappropriate for production. + ! + ! HG! + ! Added two more options to deal with forced topography changes. The topography is relaxed to a given + ! load-independent target with the chosen adjustment time scale. This may be useful to impose a topography + ! smoothly over time. A typical application would be to impose a high resolution topography after restarting from + ! an interpolated model state resulting from a lower resolution spinup. Another example would be to let a spinup + ! reach a given observed topography at present day independent of the loading history. The two options only + ! differ in how relx is initialized. Since the lithosphere is not used, it is suggested to set (lithosphere = 0). + ! To use the adjustment time scale use the default relaxing asthenosphere or set asthenosphere = 1 explicitly + ! in the [isostasy] section. 'lithosphere_period' and 'flexural_rigidity' have no impact on the results. + ! + ! - whichrelaxed = 3. Use relx as load-independent target for topg. Relax the topography to the one + ! given in relx with the chosen adjustment time scale. If relx is not specified, it is initialized to zero. + ! + ! - whichrelaxed = 4. Set relx to input topg as load-independent target for topg. Results in a time-constant + ! topography until relx is modified during the experiment. Could be used to start with the given topography + ! and then relax to a target topography that is introduced at some point as a forcing. + !------------------------------------------------------------------------- @@ -139,6 +157,48 @@ subroutine init_isostasy(model) end subroutine init_isostasy +!------------------------------------------------------------------------- + ! HG! forced topographic calculations (Greek isos "equal", zori "force") + + subroutine init_isozory(model) + + !> initialise forced topographic calculations + + use parallel + use glide_types + use glimmer_physcon, only: scyr + use glimmer_paramets, only: tim0 + implicit none + + type(glide_global_type) :: model + + model%isostasy%relaxed_tau = model%isostasy%relaxed_tau * scyr / tim0 + + end subroutine init_isozory + +!------------------------------------------------------------------------- + + subroutine isozory_compute(model) + + !> calculate forced topographic adjustment + + use glide_types + implicit none + + type(glide_global_type) :: model + + ! update bedrock abruptly if the mantle is fluid (non-viscous) + if (model%isostasy%asthenosphere == ASTHENOSPHERE_FLUID) then + model%geometry%topg = model%isostasy%relx + end if + + ! update bedrock smoothly if the mantle is relaxing + if (model%isostasy%asthenosphere == ASTHENOSPHERE_RELAXING) then + call relaxing_mantle(model) + end if + + end subroutine isozory_compute + !------------------------------------------------------------------------- subroutine isos_icewaterload(model) diff --git a/libglissade/glissade.F90 b/libglissade/glissade.F90 index 8c9f15e5..5479d6f8 100644 --- a/libglissade/glissade.F90 +++ b/libglissade/glissade.F90 @@ -244,18 +244,39 @@ subroutine glissade_initialise(model, evolve_ice) ! handle relaxed/equilibrium topo ! Initialise isostasy first - call init_isostasy(model) select case(model%isostasy%whichrelaxed) + case(RELAXED_TOPO_DEFAULT) ! topg and relx are separate input fields + + call init_isostasy(model) + ! relx and topg are used as given in input + call write_log('topg and relx are separate input fields') + case(RELAXED_TOPO_INPUT) ! supplied input topography is relaxed + call init_isostasy(model) model%isostasy%relx = model%geometry%topg + call write_log('supplied input topography is relaxed') case(RELAXED_TOPO_COMPUTE) ! supplied topography is in equilibrium !TODO - Test the case RELAXED_TOPO_COMPUTE - + call init_isostasy(model) call isos_relaxed(model) + call write_log('supplied input topography is relaxed') + + case(RELAXED_TOPO_TARGET) ! use relx as load-independent target for topg + + call init_isozory(model) + ! relx and topg are used as given in input + call write_log('use relx as load-independent target for topg') + + case(RELAXED_TOPO_FORCED) ! set relx to topg as load-independent target + ! = no isostatic adjustemnt unless relx changes + + call init_isozory(model) + model%isostasy%relx = model%geometry%topg + call write_log('set relx to topg as load-independent target') end select @@ -2129,7 +2150,11 @@ subroutine glissade_isostasy_solve(model) print*, 'Update lithospheric load: tstep_count, nlith =', & model%numerics%tstep_count, model%isostasy%nlith endif - call isos_icewaterload(model) + ! load not needed in case RELAXED_TOPO_TARGET(3) and RELAXED_TOPO_FORCED(4) + if (model%isostasy%whichrelaxed < 3) then + call isos_icewaterload(model) + end if + ! model%isostasy%new_load is used by all methods model%isostasy%new_load = .true. end if endif ! nlith > 0 @@ -2141,7 +2166,11 @@ subroutine glissade_isostasy_solve(model) ! ------------------------------------------------------------------------ if (model%options%isostasy == ISOSTASY_COMPUTE) then - call isos_compute(model) + if (model%isostasy%whichrelaxed < 3) then + call isos_compute(model) + else + call isozory_compute(model) + end if ! update topography in halo cells ! Note: For outflow BCs, most fields (thck, usrf, temp, etc.) are set to zero in the global halo,