Ability to set a diffferent reference state from the initial state should be restored

[alex-robinson]

In the latest version of FastIsostasy (5371ed), the API of isos_init_state was changed from previous versions with the set_ref argument removed. Looking at the code, I see the logic is now that either we are using a restart file, in which case the reference state is loaded from there, or the current state is assumed to represent the reference state.

However, this does not allow for the case where one would like to start from a given state that is different than the assumed reference state, but that is not obtained from a restart file. For example, this could be running an LGM NH ice sheet simulation where we start with the ICE-6G distribution, but we would still like to use present day as the reference state. This was the idea behind the previous formulation of the API.

Perhaps using the same routine isos_init_state with set_ref=.true./.false. was not the best way to handle this, and it is more challenging now with the restart capability there. Maybe the easiest is one additional public facing routine like isos_init_ref:

    subroutine isos_init_ref(isos, z_bed, H_ice, bsl, dz_ss)

        implicit none

        type(isos_class), intent(INOUT) :: isos
        real(wp), intent(IN) :: z_bed(:, :)                 ! [m] Bedrock elevation
        real(wp), intent(IN) :: H_ice(:, :)                 ! [m] Ice thickness
        real(wp), intent(IN), optional :: bsl               ! [a] Barystatic sea level
        real(wp), intent(IN), optional :: dz_ss(:, :)       ! [m] Sea surface perturbation
        
        isos%ref%z_bed(isos%domain%icrop1:isos%domain%icrop2, &
                       isos%domain%jcrop1:isos%domain%jcrop2) = z_bed
        isos%ref%Hice(isos%domain%icrop1:isos%domain%icrop2, &
                      isos%domain%jcrop1:isos%domain%jcrop2)  = H_ice
        
        if (present(bsl))   isos%ref%bsl   = bsl
        if (present(dz_ss)) isos%ref%dz_ss = dz_ss
        ! No need to set w, we, dz_ss and bsl to 0 because `zero_init_state(isos%ref)` in isos_init()
        
        ! Now the reference state has been defined
        isos%par%ref_was_set = .TRUE. 

        return

    end subroutine isos_init_ref

If this routine is not called, and no restart file is loaded, then the initial state would be used as the reference state anyway as it is now. So in isos_init_state we would have:

            if (.not. isos%par%ref_was_set) then
                call isos_init_ref(isos,isos%now%z_bed,isos%now%Hice,isos%now%bsl,isos%now%dz_ss)
                write(*,*) "isos_init_state:: reference state was set to initial state."
            end if

An alternative solution would be to add additional optional arguments to isos_init_state(..., z_bed_ref, H_ice_ref, bsl_ref, dz_ss_ref). But this seems a little bit cumbersome.