Ice-shelf solo driver and MISMIP+ updates

config_src/drivers/ice_solo_driver/ice_shelf_driver.F90

@@ -24,7 +24,7 @@ program Shelf_main
   use MOM_cpu_clock,       only : cpu_clock_id, cpu_clock_begin, cpu_clock_end
   use MOM_cpu_clock,       only : CLOCK_COMPONENT
   use MOM_debugging,       only : MOM_debugging_init
-  use MOM_diag_mediator,   only : diag_mediator_init, diag_mediator_infrastructure_init
+  use MOM_diag_mediator,   only : diag_mediator_init, diag_mediator_infrastructure_init, set_axes_info
   use MOM_diag_mediator,   only : diag_mediator_end, diag_ctrl, diag_mediator_close_registration
   use MOM_domains,         only : MOM_infra_init, MOM_infra_end
   use MOM_domains,         only : MOM_domains_init, clone_MOM_domain, pass_var
@@ -54,6 +54,8 @@ program Shelf_main
   use MOM_verticalGrid,    only : verticalGrid_type, verticalGridInit, verticalGridEnd
   use MOM_write_cputime,   only : write_cputime, MOM_write_cputime_init
   use MOM_write_cputime,   only : write_cputime_start_clock, write_cputime_CS
+  use MOM_forcing_type,    only : forcing
+  use MOM_ice_shelf_initialize, only : initialize_ice_SMB
 
   use MOM_ice_shelf, only : initialize_ice_shelf, ice_shelf_end, ice_shelf_CS
   use MOM_ice_shelf, only : ice_shelf_save_restart, solo_step_ice_shelf
@@ -75,7 +77,9 @@ program Shelf_main
   ! CPU time limit.  nmax is determined by evaluating the CPU time used between successive calls to
   ! write_cputime.  Initially it is set to be very large.
   integer :: nmax=2000000000
-
+  ! A structure containing pointers to the thermodynamic forcing fields
+  ! at the ocean surface.
+  type(forcing) :: fluxes
   ! A structure containing several relevant directory paths.
   type(directories) :: dirs
 
@@ -104,7 +108,7 @@ program Shelf_main
   real :: time_step               ! The time step [T ~> s]
 
   ! A pointer to a structure containing metrics and related information.
-  type(ocean_grid_type), pointer :: ocn_grid
+  type(ocean_grid_type), pointer :: ocn_grid => NULL()
 
   type(dyn_horgrid_type), pointer :: dG => NULL()   ! A dynamic version of the horizontal grid
   type(hor_index_type),   pointer :: HI => NULL()   ! A hor_index_type for array extents
@@ -114,7 +118,7 @@ program Shelf_main
   type(ocean_OBC_type),          pointer :: OBC => NULL()
 
   ! A pointer to a structure containing dimensional unit scaling factors.
-  type(unit_scale_type), pointer :: US
+  type(unit_scale_type), pointer :: US => NULL()
 
   type(diag_ctrl), pointer :: &
     diag => NULL()              ! A pointer to the diagnostic regulatory structure
@@ -138,8 +142,9 @@ program Shelf_main
   integer :: yr, mon, day, hr, mins, sec   ! Temp variables for writing the date.
   type(param_file_type) :: param_file      ! The structure indicating the file(s)
                                            ! containing all run-time parameters.
+  real :: smb !A constant surface mass balance that can be specified in the param_file
   character(len=9)  :: month
-  character(len=16) :: calendar = 'julian'
+  character(len=16) :: calendar = 'noleap'
   integer :: calendar_type=-1
 
   integer :: unit, io_status, ierr
@@ -184,6 +189,8 @@ program Shelf_main
     endif
   endif
 
+  ! Get the names of the I/O directories and initialization file.
+  ! Also calls the subroutine that opens run-time parameter files.
   call Get_MOM_Input(param_file, dirs)
 
   ! Read ocean_solo restart, which can override settings from the namelist.
@@ -252,8 +259,11 @@ program Shelf_main
   ! Set up the ocean model domain and grid; the ice model grid is set in initialize_ice_shelf,
   ! but the grids have strong commonalities in this configuration, and the ocean grid is required
   ! to set up the diag mediator control structure.
-  call MOM_domains_init(ocn_grid%domain, param_file)
+  allocate(ocn_grid)
+  call MOM_domains_init(ocn_grid%domain, param_file) !, domain_name='MOM')
+  allocate(HI)
   call hor_index_init(ocn_grid%Domain, HI, param_file)
+  allocate(dG)
   call create_dyn_horgrid(dG, HI)
   call clone_MOM_domain(ocn_grid%Domain, dG%Domain)
 
@@ -266,11 +276,16 @@ program Shelf_main
   ! Initialize the diag mediator.  The ocean's vertical grid is not really used here, but at
   ! present the interface to diag_mediator_init assumes the presence of ocean-specific information.
   call verticalGridInit(param_file, GV, US)
+  allocate(diag)
   call diag_mediator_init(ocn_grid, GV, US, GV%ke, param_file, diag, doc_file_dir=dirs%output_directory)
 
   call callTree_waypoint("returned from diag_mediator_init()")
 
-  call initialize_ice_shelf(param_file, ocn_grid, Time, ice_shelf_CSp, diag)
+  call set_axes_info(ocn_grid, GV, US, param_file, diag)
+
+  call initialize_ice_shelf(param_file, ocn_grid, Time, ice_shelf_CSp, diag, fluxes_in=fluxes, solo_ice_sheet_in=.true.)
+
+  call initialize_ice_SMB(fluxes%shelf_sfc_mass_flux, ocn_grid, US, param_file)
 
   ! This is the end of the code that is the counterpart of MOM_initialization.
   call callTree_waypoint("End of ice shelf initialization.")
@@ -378,7 +393,7 @@ program Shelf_main
 
     ! This call steps the model over a time time_step.
     Time1 = Master_Time ; Time = Master_Time
-    call solo_step_ice_shelf(ice_shelf_CSp, Time_step_shelf, ns_ice, Time)
+    call solo_step_ice_shelf(ice_shelf_CSp, Time_step_shelf, ns_ice, Time, fluxes_in=fluxes)
 
 !   Time = Time + Time_step_shelf
 !   This is here to enable fractional-second time steps.
@@ -412,6 +427,20 @@ program Shelf_main
       if (BTEST(Restart_control,0)) then
         call ice_shelf_save_restart(ice_shelf_CSp, Time, dirs%restart_output_dir)
       endif
+      ! Write ice shelf solo restart file.
+      if (is_root_pe())then
+        call open_ASCII_file(unit, trim(dirs%restart_output_dir)//'shelf.res')
+        write(unit, '(i6,8x,a)') calendar_type, &
+          '(Calendar: no_calendar=0, thirty_day_months=1, julian=2, gregorian=3, noleap=4)'
+
+        call get_date(Start_time, yr, mon, day, hr, mins, sec)
+        write(unit, '(6i6,8x,a)') yr, mon, day, hr, mins, sec, &
+          'Model start time:   year, month, day, hour, minute, second'
+        call get_date(Time, yr, mon, day, hr, mins, sec)
+        write(unit, '(6i6,8x,a)') yr, mon, day, hr, mins, sec, &
+          'Current model time: year, month, day, hour, minute, second'
+        call close_file(unit)
+      endif
       restart_time = restart_time + restint
     endif
 
@@ -456,12 +485,11 @@ program Shelf_main
   endif
 
   call callTree_waypoint("End Shelf_main")
+  call ice_shelf_end(ice_shelf_CSp)
   call diag_mediator_end(Time, diag, end_diag_manager=.true.)
   if (cpu_steps > 0) call write_cputime(Time, ns-1, write_CPU_CSp, call_end=.true.)
   call cpu_clock_end(termClock)
 
   call io_infra_end ; call MOM_infra_end
 
-  call ice_shelf_end(ice_shelf_CSp)
-
 end program Shelf_main

config_src/drivers/solo_driver/MOM_driver.F90

@@ -49,6 +49,7 @@ program MOM6
   use MOM_ice_shelf,       only : shelf_calc_flux, add_shelf_forces, ice_shelf_save_restart
   use MOM_ice_shelf,       only : initialize_ice_shelf_fluxes, initialize_ice_shelf_forces
   use MOM_ice_shelf,       only : ice_shelf_query
+  use MOM_ice_shelf_initialize, only : initialize_ice_SMB
   use MOM_interpolate,     only : time_interp_external_init
   use MOM_io,              only : file_exists, open_ASCII_file, close_file
   use MOM_io,              only : check_nml_error, io_infra_init, io_infra_end
@@ -134,7 +135,7 @@ program MOM6
   real :: dtdia                   ! The diabatic timestep [T ~> s]
   real :: t_elapsed_seg           ! The elapsed time in this run segment [T ~> s]
   integer :: n, ns, n_max, nts, n_last_thermo
-  logical :: diabatic_first, single_step_call
+  logical :: diabatic_first, single_step_call, initialize_smb
   type(time_type) :: Time2, time_chg ! Temporary time variables
 
   integer :: Restart_control    ! An integer that is bit-tested to determine whether
@@ -302,6 +303,9 @@ program MOM6
     call initialize_ice_shelf_forces(ice_shelf_CSp, grid, US, forces)
     call ice_shelf_query(ice_shelf_CSp, grid, data_override_shelf_fluxes=override_shelf_fluxes)
     if (override_shelf_fluxes) call data_override_init(Ocean_Domain_in=grid%domain%mpp_domain)
+    call get_param(param_file, mod_name, "INITIALIZE_ICE_SHEET_SMB", &
+                   initialize_smb, "Read in a constant SMB for the ice sheet", default=.false.)
+    if (initialize_smb) call initialize_ice_SMB(fluxes%shelf_sfc_mass_flux, grid, US, param_file)
   endif
 
 

src/framework/MOM_diag_mediator.F90

@@ -3539,37 +3539,45 @@ subroutine diag_mediator_end(time, diag_CS, end_diag_manager)
   enddo
 
   call diag_grid_storage_end(diag_cs%diag_grid_temp)
-  deallocate(diag_cs%mask3dTL)
-  deallocate(diag_cs%mask3dBL)
-  deallocate(diag_cs%mask3dCuL)
-  deallocate(diag_cs%mask3dCvL)
-  deallocate(diag_cs%mask3dTi)
-  deallocate(diag_cs%mask3dBi)
-  deallocate(diag_cs%mask3dCui)
-  deallocate(diag_cs%mask3dCvi)
+  if (associated(diag_cs%mask3dTL))  deallocate(diag_cs%mask3dTL)
+  if (associated(diag_cs%mask3dBL))  deallocate(diag_cs%mask3dBL)
+  if (associated(diag_cs%mask3dCuL)) deallocate(diag_cs%mask3dCuL)
+  if (associated(diag_cs%mask3dCvL)) deallocate(diag_cs%mask3dCvL)
+  if (associated(diag_cs%mask3dTi))  deallocate(diag_cs%mask3dTi)
+  if (associated(diag_cs%mask3dBi))  deallocate(diag_cs%mask3dBi)
+  if (associated(diag_cs%mask3dCui)) deallocate(diag_cs%mask3dCui)
+  if (associated(diag_cs%mask3dCvi)) deallocate(diag_cs%mask3dCvi)
   do dl=2,MAX_DSAMP_LEV
-    deallocate(diag_cs%dsamp(dl)%mask2dT)
-    deallocate(diag_cs%dsamp(dl)%mask2dBu)
-    deallocate(diag_cs%dsamp(dl)%mask2dCu)
-    deallocate(diag_cs%dsamp(dl)%mask2dCv)
-    deallocate(diag_cs%dsamp(dl)%mask3dTL)
-    deallocate(diag_cs%dsamp(dl)%mask3dBL)
-    deallocate(diag_cs%dsamp(dl)%mask3dCuL)
-    deallocate(diag_cs%dsamp(dl)%mask3dCvL)
-    deallocate(diag_cs%dsamp(dl)%mask3dTi)
-    deallocate(diag_cs%dsamp(dl)%mask3dBi)
-    deallocate(diag_cs%dsamp(dl)%mask3dCui)
-    deallocate(diag_cs%dsamp(dl)%mask3dCvi)
+    if (associated(diag_cs%dsamp(dl)%mask2dT))   deallocate(diag_cs%dsamp(dl)%mask2dT)
+    if (associated(diag_cs%dsamp(dl)%mask2dBu))  deallocate(diag_cs%dsamp(dl)%mask2dBu)
+    if (associated(diag_cs%dsamp(dl)%mask2dCu))  deallocate(diag_cs%dsamp(dl)%mask2dCu)
+    if (associated(diag_cs%dsamp(dl)%mask2dCv))  deallocate(diag_cs%dsamp(dl)%mask2dCv)
+    if (associated(diag_cs%dsamp(dl)%mask3dTL))  deallocate(diag_cs%dsamp(dl)%mask3dTL)
+    if (associated(diag_cs%dsamp(dl)%mask3dBL))  deallocate(diag_cs%dsamp(dl)%mask3dBL)
+    if (associated(diag_cs%dsamp(dl)%mask3dCuL)) deallocate(diag_cs%dsamp(dl)%mask3dCuL)
+    if (associated(diag_cs%dsamp(dl)%mask3dCvL)) deallocate(diag_cs%dsamp(dl)%mask3dCvL)
+    if (associated(diag_cs%dsamp(dl)%mask3dTi))  deallocate(diag_cs%dsamp(dl)%mask3dTi)
+    if (associated(diag_cs%dsamp(dl)%mask3dBi))  deallocate(diag_cs%dsamp(dl)%mask3dBi)
+    if (associated(diag_cs%dsamp(dl)%mask3dCui)) deallocate(diag_cs%dsamp(dl)%mask3dCui)
+    if (associated(diag_cs%dsamp(dl)%mask3dCvi)) deallocate(diag_cs%dsamp(dl)%mask3dCvi)
 
     do i=1,diag_cs%num_diag_coords
-      deallocate(diag_cs%dsamp(dl)%remap_axesTL(i)%dsamp(dl)%mask3d)
-      deallocate(diag_cs%dsamp(dl)%remap_axesCuL(i)%dsamp(dl)%mask3d)
-      deallocate(diag_cs%dsamp(dl)%remap_axesCvL(i)%dsamp(dl)%mask3d)
-      deallocate(diag_cs%dsamp(dl)%remap_axesBL(i)%dsamp(dl)%mask3d)
-      deallocate(diag_cs%dsamp(dl)%remap_axesTi(i)%dsamp(dl)%mask3d)
-      deallocate(diag_cs%dsamp(dl)%remap_axesCui(i)%dsamp(dl)%mask3d)
-      deallocate(diag_cs%dsamp(dl)%remap_axesCvi(i)%dsamp(dl)%mask3d)
-      deallocate(diag_cs%dsamp(dl)%remap_axesBi(i)%dsamp(dl)%mask3d)
+      if (associated(diag_cs%dsamp(dl)%remap_axesTL(i)%dsamp(dl)%mask3d)) &
+        deallocate(diag_cs%dsamp(dl)%remap_axesTL(i)%dsamp(dl)%mask3d)
+      if (associated(diag_cs%dsamp(dl)%remap_axesCuL(i)%dsamp(dl)%mask3d)) &
+        deallocate(diag_cs%dsamp(dl)%remap_axesCuL(i)%dsamp(dl)%mask3d)
+      if (associated(diag_cs%dsamp(dl)%remap_axesCvL(i)%dsamp(dl)%mask3d)) &
+        deallocate(diag_cs%dsamp(dl)%remap_axesCvL(i)%dsamp(dl)%mask3d)
+      if (associated(diag_cs%dsamp(dl)%remap_axesBL(i)%dsamp(dl)%mask3d)) &
+        deallocate(diag_cs%dsamp(dl)%remap_axesBL(i)%dsamp(dl)%mask3d)
+      if (associated(diag_cs%dsamp(dl)%remap_axesTi(i)%dsamp(dl)%mask3d)) &
+        deallocate(diag_cs%dsamp(dl)%remap_axesTi(i)%dsamp(dl)%mask3d)
+      if (associated(diag_cs%dsamp(dl)%remap_axesCui(i)%dsamp(dl)%mask3d)) &
+        deallocate(diag_cs%dsamp(dl)%remap_axesCui(i)%dsamp(dl)%mask3d)
+      if (associated(diag_cs%dsamp(dl)%remap_axesCvi(i)%dsamp(dl)%mask3d)) &
+        deallocate(diag_cs%dsamp(dl)%remap_axesCvi(i)%dsamp(dl)%mask3d)
+      if (associated(diag_cs%dsamp(dl)%remap_axesBi(i)%dsamp(dl)%mask3d)) &
+        deallocate(diag_cs%dsamp(dl)%remap_axesBi(i)%dsamp(dl)%mask3d)
     enddo
   enddo
 

src/ice_shelf/MOM_ice_shelf.F90

@@ -1873,7 +1873,8 @@ subroutine initialize_ice_shelf_fluxes(CS, ocn_grid, US, fluxes_in)
          tau_mag=.true.)
   else
     call MOM_mesg("MOM_ice_shelf.F90, initialize_ice_shelf: allocating fluxes in solo mode.")
-    call allocate_forcing_type(CS%Grid_in, fluxes_in, ustar=.true., shelf=.true., press=.true., tau_mag=.true.)
+    call allocate_forcing_type(CS%Grid_in, fluxes_in, ustar=.true., shelf=.true., &
+         press=.true., shelf_sfc_accumulation = CS%active_shelf_dynamics, tau_mag=.true.)
   endif
   if (CS%rotate_index) then
     allocate(fluxes)
@@ -2178,20 +2179,22 @@ subroutine ice_shelf_end(CS)
 end subroutine ice_shelf_end
 
 !> This routine is for stepping a stand-alone ice shelf model without an ocean.
-subroutine solo_step_ice_shelf(CS, time_interval, nsteps, Time, min_time_step_in)
+subroutine solo_step_ice_shelf(CS, time_interval, nsteps, Time, min_time_step_in, fluxes_in)
   type(ice_shelf_CS), pointer    :: CS      !< A pointer to the ice shelf control structure
   type(time_type), intent(in)    :: time_interval !< The time interval for this update [s].
   integer,         intent(inout) :: nsteps  !< The running number of ice shelf steps.
   type(time_type), intent(inout) :: Time    !< The current model time
   real,  optional, intent(in)    :: min_time_step_in !< The minimum permitted time step [T ~> s].
-
+  type(forcing),      optional, target, intent(inout) :: fluxes_in !< A structure containing pointers to any
+                                                         !!  possible thermodynamic or mass-flux forcing fields.
   type(ocean_grid_type), pointer :: G => NULL()  ! A pointer to the ocean's grid structure
   type(unit_scale_type), pointer :: US => NULL() ! Pointer to a structure containing
                                                  ! various unit conversion factors
   type(ice_shelf_state), pointer :: ISS => NULL() !< A structure with elements that describe
                                           !! the ice-shelf state
   real :: remaining_time    ! The remaining time in this call [T ~> s]
   real :: time_step         ! The internal time step during this call [T ~> s]
+  real :: full_time_step    ! The external time step (sum of internal time steps) during this call [T ~> s]
   real :: min_time_step     ! The minimal required timestep that would indicate a fatal problem [T ~> s]
   character(len=240) :: mesg
   logical :: update_ice_vel ! If true, it is time to update the ice shelf velocities.
@@ -2205,6 +2208,7 @@ subroutine solo_step_ice_shelf(CS, time_interval, nsteps, Time, min_time_step_in
   is = G%isc ; iec = G%iec ; js = G%jsc ; jec = G%jec
 
   remaining_time = US%s_to_T*time_type_to_real(time_interval)
+  full_time_step = remaining_time
 
   if (present (min_time_step_in)) then
     min_time_step = min_time_step_in
@@ -2228,6 +2232,8 @@ subroutine solo_step_ice_shelf(CS, time_interval, nsteps, Time, min_time_step_in
       call MOM_mesg("solo_step_ice_shelf: "//mesg, 5)
     endif
 
+    call change_thickness_using_precip(CS, ISS, G, US, fluxes_in, time_step, Time)
+
     remaining_time = remaining_time - time_step
 
     ! If the last mini-timestep is a day or less, we cannot expect velocities to change by much.
@@ -2237,13 +2243,13 @@ subroutine solo_step_ice_shelf(CS, time_interval, nsteps, Time, min_time_step_in
 
     call update_ice_shelf(CS%dCS, ISS, G, US, time_step, Time, must_update_vel=update_ice_vel)
 
-    call enable_averages(time_step, Time, CS%diag)
+  enddo
+
+  call enable_averages(full_time_step, Time, CS%diag)
     if (CS%id_area_shelf_h > 0) call post_data(CS%id_area_shelf_h, ISS%area_shelf_h, CS%diag)
     if (CS%id_h_shelf > 0) call post_data(CS%id_h_shelf, ISS%h_shelf, CS%diag)
     if (CS%id_h_mask > 0) call post_data(CS%id_h_mask, ISS%hmask, CS%diag)
-    call disable_averaging(CS%diag)
-
-  enddo
+  call disable_averaging(CS%diag)
 
 end subroutine solo_step_ice_shelf
 

src/ice_shelf/MOM_ice_shelf_dynamics.F90

@@ -768,7 +768,7 @@ subroutine update_ice_shelf(CS, ISS, G, US, time_step, Time, ocean_mass, coupled
 
 ! call ice_shelf_temp(CS, ISS, G, US, time_step, ISS%water_flux, Time)
 
-  if (update_ice_vel) then
+  if (CS%elapsed_velocity_time >= CS%velocity_update_time_step) then
     call enable_averages(CS%elapsed_velocity_time, Time, CS%diag)
     if (CS%id_col_thick > 0) call post_data(CS%id_col_thick, CS%OD_av, CS%diag)
     if (CS%id_u_shelf > 0) call post_data(CS%id_u_shelf, CS%u_shelf, CS%diag)

src/ice_shelf/MOM_ice_shelf_initialize.F90

@@ -22,6 +22,7 @@ module MOM_ice_shelf_initialize
 public initialize_ice_shelf_boundary_from_file
 public initialize_ice_C_basal_friction
 public initialize_ice_AGlen
+public initialize_ice_SMB
 ! A note on unit descriptions in comments: MOM6 uses units that can be rescaled for dimensional
 ! consistency testing. These are noted in comments with units like Z, H, L, and T, along with
 ! their mks counterparts with notation like "a velocity [Z T-1 ~> m s-1]".  If the units
@@ -657,5 +658,51 @@ subroutine initialize_ice_AGlen(AGlen, G, US, PF)
     call MOM_read_data(filename,trim(varname), AGlen, G%Domain)
 
   endif
-end subroutine
+end subroutine initialize_ice_AGlen
+
+!> Initialize ice surface mass balance field that is held constant over time
+subroutine initialize_ice_SMB(SMB, G, US, PF)
+  type(ocean_grid_type), intent(in)    :: G    !< The ocean's grid structure
+  real, dimension(SZDI_(G),SZDJ_(G)), &
+                         intent(inout) :: SMB !< Ice surface mass balance parameter, often in [kg m-2 s-1]
+  type(unit_scale_type), intent(in)    :: US !< A structure containing unit conversion factors
+  type(param_file_type), intent(in)    :: PF !< A structure to parse for run-time parameters
+
+  real :: SMB_val  ! Constant ice surface mass balance parameter, often in [kg m-2 s-1]
+  character(len=40)  :: mdl = "initialize_ice_SMB" ! This subroutine's name.
+  character(len=200) :: config
+  character(len=200) :: varname
+  character(len=200) :: inputdir, filename, SMB_file
+
+  call get_param(PF, mdl, "ICE_SMB_CONFIG", config, &
+                 "This specifies how the initial ice surface mass balance parameter is specified. "//&
+                 "Valid values are: CONSTANT and FILE.", &
+                 default="CONSTANT")
+
+  if (trim(config)=="CONSTANT") then
+    call get_param(PF, mdl, "SMB", SMB_val, &
+                 "Surface mass balance.", units="kg m-2 s-1", default=0.0)
+
+    SMB(:,:) = SMB_val
+
+  elseif (trim(config)=="FILE") then
+    call MOM_mesg("  MOM_ice_shelf.F90, initialize_ice_shelf: reading SMB parameter")
+    call get_param(PF, mdl, "INPUTDIR", inputdir, default=".")
+    inputdir = slasher(inputdir)
+
+    call get_param(PF, mdl, "ICE_SMB_FILE", SMB_file, &
+                 "The file from which the ice surface mass balance is read.", &
+                 default="ice_SMB.nc")
+    filename = trim(inputdir)//trim(SMB_file)
+    call log_param(PF, mdl, "INPUTDIR/ICE_SMB_FILE", filename)
+    call get_param(PF, mdl, "ICE_SMB_VARNAME", varname, &
+                   "The variable to use as surface mass balance.", &
+                   default="SMB")
+
+    if (.not.file_exists(filename, G%Domain)) call MOM_error(FATAL, &
+       " initialize_ice_SMV_from_file: Unable to open "//trim(filename))
+    call MOM_read_data(filename,trim(varname), SMB, G%Domain)
+
+  endif
+end subroutine initialize_ice_SMB
 end module MOM_ice_shelf_initialize