Adding forced topography options

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

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

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)

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,