Changes to fix merge inconsistancies after rebase.

Author: Katetc

.gitignore

@@ -43,3 +43,25 @@ P75R
 P75S
 Stnd
 mismip+Ice0
+
+#Ignore compiled topography files
+bin_to_cube/bin_to_cube
+bin_to_cube/bin_to_cube.o
+bin_to_cube/shr_kind_mod.mod
+bin_to_cube/shr_kind_mod.o
+cube_to_target/cube_to_target
+cube_to_target/cube_to_target.o
+cube_to_target/reconstruct.mod
+cube_to_target/reconstruct.o
+cube_to_target/remap.mod
+cube_to_target/remap.o
+cube_to_target/ridge_ana.mod
+cube_to_target/ridge_ana.o
+cube_to_target/rot.o
+cube_to_target/rotation.mod
+cube_to_target/shared_vars.mod
+cube_to_target/shared_vars.o
+cube_to_target/shr_kind_mod.mod
+cube_to_target/shr_kind_mod.o
+cube_to_target/smooth_topo_cube.o
+cube_to_target/smooth_topo_cube_sph.mod

libglide/glide_setup.F90

@@ -764,6 +764,7 @@ subroutine handle_options(section, model)
     call GetValue(section,'forcewrite_final', model%options%forcewrite_final)
     call GetValue(section,'restart',model%options%is_restart)
     call GetValue(section,'restart_extend_velo',model%options%restart_extend_velo)
+    call GetValue(section,'forcewrite_restart',model%options%forcewrite_restart)
 
   end subroutine handle_options
 
@@ -1666,6 +1667,17 @@ subroutine print_options(model)
        endif
     end if
 
+    if (model%options%forcewrite_restart) then
+       call write_log('Will write to output files on restart')
+    endif
+
+!!     This option is not currently supported
+!!    if (model%options%which_bproc < 0 .or. model%options%which_bproc >= size(which_bproc)) then
+!!       call write_log('Error, basal_proc out of range',GM_FATAL)
+!!    end if
+!!    write(message,*) 'basal_proc              : ',model%options%which_bproc,which_bproc(model%options%which_bproc)
+!!    call write_log(message)
+
     !HO options
 
     if (model%options%whichdycore /= DYCORE_GLIDE) then   ! glissade higher-order
@@ -2640,10 +2652,6 @@ subroutine print_parameters(model)
     elseif (model%options%which_ho_babc == HO_BABC_POWERLAW) then
        write(message,*) 'Cp for power law, Pa (m/yr)^(-1/3)           : ', model%basal_physics%powerlaw_c_const
        call write_log(message)
-       write(message,*) 'Max Cp for power law, Pa (m/yr)^(-1/3)       : ', model%basal_physics%powerlaw_c_max
-       call write_log(message)
-       write(message,*) 'Min Cp for power law, Pa (m/yr)^(-1/3)       : ', model%basal_physics%powerlaw_c_min
-       call write_log(message)
        write(message,*) 'm exponent for power law                     : ', model%basal_physics%powerlaw_m
        call write_log(message)
     elseif (model%options%which_ho_babc == HO_BABC_COULOMB_FRICTION) then
@@ -2660,10 +2668,6 @@ subroutine print_parameters(model)
        call write_log(message)
        write(message,*) 'Cp for Schoof power law, Pa (m/yr)^(-1/3)    : ', model%basal_physics%powerlaw_c_const
        call write_log(message)
-       write(message,*) 'Max Cp for power law, Pa (m/yr)^(-1/3)       : ', model%basal_physics%powerlaw_c_max
-       call write_log(message)
-       write(message,*) 'Min Cp for power law, Pa (m/yr)^(-1/3)       : ', model%basal_physics%powerlaw_c_min
-       call write_log(message)
        write(message,*) 'm exponent for Schoof power law              : ', model%basal_physics%powerlaw_m
        call write_log(message)
     elseif (model%options%which_ho_babc == HO_BABC_COULOMB_POWERLAW_TSAI) then
@@ -2673,10 +2677,6 @@ subroutine print_parameters(model)
        call write_log(message)
        write(message,*) 'Cp for Tsai power law, Pa (m/yr)^(-1/3)      : ', model%basal_physics%powerlaw_c_const
        call write_log(message)
-       write(message,*) 'Max Cp for power law, Pa (m/yr)^(-1/3)       : ', model%basal_physics%powerlaw_c_max
-       call write_log(message)
-       write(message,*) 'Min Cp for power law, Pa (m/yr)^(-1/3)       : ', model%basal_physics%powerlaw_c_min
-       call write_log(message)
        write(message,*) 'm exponent for Tsai power law                : ', model%basal_physics%powerlaw_m
        call write_log(message)
     elseif (model%options%which_ho_babc == HO_BABC_POWERLAW_EFFECPRESS) then
@@ -3399,6 +3399,7 @@ subroutine define_glide_restart_variables(model, model_id)
     ! Subroutine arguments
     !------------------------------------------------------------------------------------
     type(glide_global_type), intent (in) :: model  !> Derived type holding all model info
+
     integer, intent(in) :: model_id  !> identifier of this ice sheet instance (1 - N, where N is the total number of ice sheet models in this run)
 
     !------------------------------------------------------------------------------------
@@ -3458,7 +3459,7 @@ subroutine define_glide_restart_variables(model, model_id)
              call glide_add_to_restart_variable_list('smb_gradz', model_id)
           end select
 
-          call glide_add_to_restart_variable_list('smb_reference_usrf', model_id)
+          call glide_add_to_restart_variable_list('usrf_ref', model_id)
 
        case(SMB_INPUT_FUNCTION_XYZ)
 
@@ -3484,7 +3485,7 @@ subroutine define_glide_restart_variables(model, model_id)
           if (options%smb_input_function == SMB_INPUT_FUNCTION_XY_GRADZ) then
              ! usrf_ref was added to restart above; nothing to do here
           else
-             call glide_add_to_restart_variable_list('smb_reference_usrf', model_id)
+             call glide_add_to_restart_variable_list('usrf_ref', model_id)
           endif
 
        case(ARTM_INPUT_FUNCTION_XYZ)
@@ -3811,6 +3812,36 @@ subroutine define_glide_restart_variables(model, model_id)
        call glide_add_to_restart_variable_list('f_effecpress_ocean_p', model_id)
     endif
 
+    ! fields needed for inversion options that try to match local thickness or upper surface elevation
+    ! Note: If usrf_obs is supplied, thck_obs will be computed at initialization
+    if (options%which_ho_powerlaw_c == HO_POWERLAW_C_INVERSION .or. &
+        options%which_ho_coulomb_c  == HO_COULOMB_C_INVERSION  .or. &
+        options%which_ho_deltaT_ocn == HO_DELTAT_OCN_INVERSION) then
+       call glide_add_to_restart_variable_list('usrf_obs', model_id)
+       !WHL - velo_sfc_obs is not strictly needed unless inverting for surface velo,
+       !      but is handy for diagnostics
+       call glide_add_to_restart_variable_list('velo_sfc_obs', model_id)
+    endif
+
+    ! fields needed for inversion options that try to match local dthck_dt
+    ! Note: This is not strictly needed for all options, but still is a useful diagnostic.
+    if (options%which_ho_deltaT_ocn /= HO_DELTAT_OCN_NONE) then
+       call glide_add_to_restart_variable_list('dthck_dt_obs', model_id)
+       call glide_add_to_restart_variable_list('dthck_dt_obs_basin', model_id)
+    endif
+
+    ! effective pressure options
+    ! f_effecpress_bwat represents the reduction of overburden pressure from bwatflx
+    if (options%which_ho_effecpress == HO_EFFECPRESS_BWATFLX) then
+       call glide_add_to_restart_variable_list('f_effecpress_bwat', model_id)
+    endif
+
+    ! f_effecpress_ocean_p represents the reduction of overburden pressure when ocean_p > 0
+    ! Needs to be saved in case this fraction is relaxed over time toward (1 - Hf/H)^p
+    if (model%basal_physics%p_ocean_penetration > 0.0d0) then
+       call glide_add_to_restart_variable_list('f_effecpress_ocean_p', model_id)
+    endif
+
     ! geothermal heat flux option
     select case (options%gthf)
       case(GTHF_COMPUTE)
@@ -3876,6 +3907,11 @@ subroutine define_glide_restart_variables(model, model_id)
        call glide_add_to_restart_variable_list('glacier_area_init', model_id)
     endif
 
+    ! basal processes module - requires tauf for a restart
+!!    if (options%which_bproc /= BAS_PROC_DISABLED ) then
+!!        call glide_add_to_restart_variable_list('tauf', model_id)
+!!    endif
+
     ! TODO bmlt was set as a restart variable, but I'm not sure when or if it is needed.
 
     ! TODO age should be a restart variable if it is an input variable.  

libglide/glide_types.F90

@@ -774,6 +774,9 @@ module glide_types
     !>          (required if restart velocities are nonzero on global boundaries)
     !> \end{description}
 
+    logical :: forcewrite_restart = .false.
+    !> flag that indicates whether to force writing of output on restart
+
     ! This is a Glimmer serial option
     ! The parallel code enforces periodic EW and NS boundary conditions by default
     logical :: periodic_ew = .false.

libglissade/glissade.F90

@@ -398,6 +398,31 @@ subroutine glissade_initialise(model, evolve_ice)
        call check_fill_values(model%geometry%dthck_dt_obs)
     endif
 
+    ! Some input fields may have a netCDF fill value, typically a very large positive number.
+    ! If present, convert these values to zero (or optionally, another suitable value).
+    ! Note: Optionally, can pass a user-specified fill value and replacement value,
+    !        and return a mask of grid cells where values are replaced.
+    !       Depending on the input dataset, might have fill values in other fields (e.g., artm, topg)
+
+    if (model%options%smb_input == SMB_INPUT_MMYR_WE) then
+       call check_fill_values(model%climate%smb)
+    else
+       call check_fill_values(model%climate%acab)
+    endif
+
+    if (model%options%gthf == GTHF_PRESCRIBED_2D) then
+       call check_fill_values(model%temper%bheatflx)
+    endif
+
+    if (associated(model%ocean_data%thermal_forcing)) then
+       call check_fill_values(model%ocean_data%thermal_forcing)
+    endif
+
+    if (model%options%which_ho_deltaT_ocn == HO_DELTAT_OCN_DTHCK_DT .or.  &
+        model%options%enable_acab_dthck_dt_correction) then
+       call check_fill_values(model%geometry%dthck_dt_obs)
+    endif
+
     ! Allocate mask arrays in case they are needed below
     allocate(ice_mask(model%general%ewn, model%general%nsn))
     allocate(floating_mask(model%general%ewn, model%general%nsn))
@@ -936,6 +961,7 @@ subroutine glissade_initialise(model, evolve_ice)
     ! Initialize the effective pressure calculation
 
     if (model%options%is_restart == NO_RESTART) then
+
        call glissade_init_effective_pressure(model%options%which_ho_effecpress,  &
                                              model%basal_physics)
     endif
@@ -2516,7 +2542,7 @@ subroutine glissade_thickness_tracer_solve(model)
           !        a suite of automated stability tests, e.g. with the stabilitySlab.py script.
           if (advective_cfl > 1.0d0) then
              if (main_task) print*, 'advective CFL violation; call glide_finalise and exit cleanly'
-             call glide_finalise(model)
+             call glide_finalise(model, forcewrite_arg=.true.)
              stop
           else
              nsubcyc = model%numerics%subcyc
@@ -4099,11 +4125,7 @@ subroutine glissade_diagnostic_variable_solve(model)
          f_ground_cell_obs,  & ! f_ground_cell as a function of thck_obs (instead of current thck)
          f_ground_obs,       & ! f_ground as a function of thck_obs (instead of current thck)
          f_flotation_obs,    & ! f_flotation_obs as a function of thck_obs (instead of current thck)
-         thck_calving_front,   & ! effective thickness of ice at the calving front
-         powerlaw_c_icegrid      ! powerlaw_c on the unstaggered ice grid
-
-    real(dp), dimension(model%general%ewn, model%general%nsn) ::  &
-         flow_enhancement_factor_float    ! flow enhancement factor for floating ice
+         thck_calving_front    ! effective thickness of ice at the calving front
 
     real(dp) :: &
          dsigma,                   & ! layer thickness in sigma coordinates
@@ -4327,6 +4349,7 @@ subroutine glissade_diagnostic_variable_solve(model)
 
     ! Compute the thickness tendency dH/dt from one step to the next (m/s)
     ! This tendency is used for coulomb_c and powerlaw_c inversion.
+
     if ( (model%options%is_restart == STANDARD_RESTART .or. model%options%is_restart == HYBRID_RESTART) &
          .and. (model%numerics%time == model%numerics%tstart) ) then
        ! first call after a restart; do not compute dthck_dt
@@ -4528,7 +4551,6 @@ subroutine glissade_diagnostic_variable_solve(model)
 
     endif   ! which_ho_flow_enhancement_factor
 
-
     ! If glaciers are enabled, then do various updates:
     ! (1) If inverting for mu_star, alpha_snow, or powerlaw_c, then
     !     (a) Accumulate the fields needed for the inversion.

libglissade/glissade_basal_traction.F90

@@ -145,6 +145,7 @@ subroutine calcbeta (whichbabc,                    &
 
   ! variables for Coulomb friction law
   real(dp) :: coulomb_c   ! Coulomb law friction coefficient (unitless)
+  real(dp) :: powerlaw_c_const  ! power law friction coefficient (Pa m^{-1/3} yr^{1/3})
   real(dp) :: lambda_max        ! wavelength of bedrock bumps at subgrid scale (m)
   real(dp) :: m_max             ! maximum bed obstacle slope (unitless)
   real(dp) :: m                 ! exponent m in power law
@@ -210,6 +211,7 @@ subroutine calcbeta (whichbabc,                    &
                                   basal_physics%coulomb_c_bedmin,  &
                                   basal_physics%coulomb_c_bedmax,  &
                                   basal_physics%coulomb_c)
+
   endif
 
   ! Compute beta based on whichbabc

libglissade/glissade_velo_higher.F90

@@ -764,6 +764,7 @@ subroutine glissade_velo_higher_solve(model,                &
     integer ::   &
        whichbabc, &             ! option for basal boundary condition
        whichbeta_limit, &       ! option to limit beta for grounded ice
+       which_powerlaw_c, &      ! option for powerlaw friction parameter Cp
        which_coulomb_c, &       ! option for coulomb friction parameter Cc
        whichefvs, &             ! option for effective viscosity calculation 
                                 ! (calculate it or make it uniform)