Ice-shelf bugfixes for restarts and halo updates

This commit fixes a combination of bugs related to halo update errors, improperly-defined iteration bounds, and an unexpected change in ice-shelf bed_elev after a restart. These issues caused crashes during the ice-shelf velocity (SSA) solution, and ice shelf dynamics restarts were not bitwise identical. The issues were resolved with the following changes: -Added bed_elev to the ice shelf restart file (previously, bed_elev was instead calculated incorrectly after a restart). Also fixed the subsequent bed_elev pass_var call, where the halo update was failing due to an error in how an optional argument (CENTER) was passed. -Added additional pass_var/pass_vector calls upon ice shelf restarts to guarantee that all halos are properly filled. -Disabled the additional ice shelf flow solve (ice_shelf_solve_outer) after a restart, which was unnecessary and prevented bitwise identical restarts. -Fixed incorrect bounds for iterations within ice_shelf_solve_inner and for the basis functions in calc_shelf_visc This commit also sets the mol_wt argument to optional in the dummy function aof_set_coupler_flux (from ice_solo_driver/atmost_ocean_fluxes.F90), as that is how it is called from the FMS coupler. This is required for compilation in ice-shelf-only mode.
NOAA-GFDL · Sep 6, 2023 · 1577ae1 · 1577ae1
1 parent 9f7f86d
commit 1577ae1
Show file tree

Hide file tree

Showing 3 changed files with 48 additions and 33 deletions.
diff --git a/config_src/drivers/ice_solo_driver/atmos_ocean_fluxes.F90 b/config_src/drivers/ice_solo_driver/atmos_ocean_fluxes.F90
@@ -13,7 +13,7 @@ module atmos_ocean_fluxes_mod
 !> This subroutine duplicates an interface used by the FMS coupler, but only
 !! returns a value of -1.  None of the arguments are used for anything.
 function aof_set_coupler_flux(name, flux_type, implementation, atm_tr_index,     &
-                              param, flag, ice_restart_file, ocean_restart_file, &
+                              param, flag, mol_wt, ice_restart_file, ocean_restart_file, &
                               units, caller, verbosity)  result (coupler_index)
 
   character(len=*),                intent(in) :: name !< An unused argument
@@ -22,6 +22,7 @@ function aof_set_coupler_flux(name, flux_type, implementation, atm_tr_index,
   integer,               optional, intent(in) :: atm_tr_index !< An unused argument
   real,    dimension(:), optional, intent(in) :: param !< An unused argument
   logical, dimension(:), optional, intent(in) :: flag !< An unused argument
+  real,                  optional, intent(in) :: mol_wt !< An unused argument
   character(len=*),      optional, intent(in) :: ice_restart_file !< An unused argument
   character(len=*),      optional, intent(in) :: ocean_restart_file !< An unused argument
   character(len=*),      optional, intent(in) :: units !< An unused argument

diff --git a/src/ice_shelf/MOM_ice_shelf_dynamics.F90 b/src/ice_shelf/MOM_ice_shelf_dynamics.F90
@@ -265,23 +265,23 @@ subroutine register_ice_shelf_dyn_restarts(G, US, param_file, CS, restart_CS)
     call get_param(param_file, mdl, "MISSING_SHELF_TEMPERATURE", T_shelf_missing, &
                  "An ice shelf temperature to use where there is no ice shelf.",&
                  units="degC", default=-10.0, scale=US%degC_to_C, do_not_log=.true.)
-    allocate( CS%u_shelf(IsdB:IedB,JsdB:JedB), source=0.0 )
-    allocate( CS%v_shelf(IsdB:IedB,JsdB:JedB), source=0.0 )
-    allocate( CS%t_shelf(isd:ied,jsd:jed), source=T_shelf_missing ) ! [C ~> degC]
-    allocate( CS%ice_visc(isd:ied,jsd:jed), source=0.0 )
-    allocate( CS%AGlen_visc(isd:ied,jsd:jed), source=2.261e-25 ) ! [Pa-3 s-1]
-    allocate( CS%basal_traction(isd:ied,jsd:jed), source=0.0 )   ! [R L2 T-2 ~> Pa]
-    allocate( CS%C_basal_friction(isd:ied,jsd:jed), source=5.0e10 ) ! [Pa (m-1 s)^n_sliding]
-    allocate( CS%OD_av(isd:ied,jsd:jed), source=0.0 )
-    allocate( CS%ground_frac(isd:ied,jsd:jed), source=0.0 )
-    allocate( CS%taudx_shelf(IsdB:IedB,JsdB:JedB), source=0.0 )
-    allocate( CS%taudy_shelf(IsdB:IedB,JsdB:JedB), source=0.0 )
-    allocate( CS%bed_elev(isd:ied,jsd:jed) ) ; CS%bed_elev(:,:) = G%bathyT(:,:) + G%Z_ref
-    allocate( CS%u_bdry_val(IsdB:IedB,JsdB:JedB), source=0.0 )
-    allocate( CS%v_bdry_val(IsdB:IedB,JsdB:JedB), source=0.0 )
-    allocate( CS%u_face_mask_bdry(IsdB:IedB,JsdB:JedB), source=-2.0 )
-    allocate( CS%v_face_mask_bdry(IsdB:iedB,JsdB:JedB), source=-2.0 )
-    allocate( CS%h_bdry_val(isd:ied,jsd:jed), source=0.0 )
+    allocate(CS%u_shelf(IsdB:IedB,JsdB:JedB), source=0.0)
+    allocate(CS%v_shelf(IsdB:IedB,JsdB:JedB), source=0.0)
+    allocate(CS%t_shelf(isd:ied,jsd:jed), source=T_shelf_missing) ! [C ~> degC]
+    allocate(CS%ice_visc(isd:ied,jsd:jed), source=0.0)
+    allocate(CS%AGlen_visc(isd:ied,jsd:jed), source=2.261e-25) ! [Pa-3 s-1]
+    allocate(CS%basal_traction(isd:ied,jsd:jed), source=0.0)   ! [R L2 T-2 ~> Pa]
+    allocate(CS%C_basal_friction(isd:ied,jsd:jed), source=5.0e10) ! [Pa (m-1 s)^n_sliding]
+    allocate(CS%OD_av(isd:ied,jsd:jed), source=0.0)
+    allocate(CS%ground_frac(isd:ied,jsd:jed), source=0.0)
+    allocate(CS%taudx_shelf(IsdB:IedB,JsdB:JedB), source=0.0)
+    allocate(CS%taudy_shelf(IsdB:IedB,JsdB:JedB), source=0.0)
+    allocate(CS%bed_elev(isd:ied,jsd:jed), source=0.0)
+    allocate(CS%u_bdry_val(IsdB:IedB,JsdB:JedB), source=0.0)
+    allocate(CS%v_bdry_val(IsdB:IedB,JsdB:JedB), source=0.0)
+    allocate(CS%u_face_mask_bdry(IsdB:IedB,JsdB:JedB), source=-2.0)
+    allocate(CS%v_face_mask_bdry(IsdB:iedB,JsdB:JedB), source=-2.0)
+    allocate(CS%h_bdry_val(isd:ied,jsd:jed), source=0.0)
    ! additional restarts for ice shelf state
     call register_restart_field(CS%u_shelf, "u_shelf", .false., restart_CS, &
                                 "ice sheet/shelf u-velocity", &
@@ -310,6 +310,8 @@ subroutine register_ice_shelf_dyn_restarts(G, US, param_file, CS, restart_CS)
                                 "ice-stiffness parameter", "Pa-3 s-1")
     call register_restart_field(CS%h_bdry_val, "h_bdry_val", .false., restart_CS, &
                                 "ice thickness at the boundary", "m", conversion=US%Z_to_m)
+    call register_restart_field(CS%bed_elev, "bed elevation", .true., restart_CS, &
+                                "bed elevation", "m", conversion=US%Z_to_m)
   endif
 
 end subroutine register_ice_shelf_dyn_restarts
@@ -509,7 +511,15 @@ subroutine initialize_ice_shelf_dyn(param_file, Time, ISS, CS, G, US, diag, new_
     call pass_var(CS%ice_visc,G%domain)
     call pass_var(CS%basal_traction, G%domain)
     call pass_var(CS%AGlen_visc, G%domain)
+    call pass_var(CS%bed_elev, G%domain)
+    call pass_var(CS%C_basal_friction, G%domain)
+    call pass_var(CS%h_bdry_val, G%domain)
+    call pass_var(CS%thickness_bdry_val, G%domain)
+
     call pass_vector(CS%u_shelf, CS%v_shelf, G%domain, TO_ALL, BGRID_NE)
+    call pass_vector(CS%u_bdry_val, CS%v_bdry_val, G%domain, TO_ALL, BGRID_NE)
+    call pass_vector(CS%u_face_mask_bdry, CS%v_face_mask_bdry, G%domain, TO_ALL, BGRID_NE)
+    call update_velocity_masks(CS, G, ISS%hmask, CS%umask, CS%vmask, CS%u_face_mask, CS%v_face_mask)
   endif
 
   if (active_shelf_dynamics) then
@@ -561,7 +571,8 @@ subroutine initialize_ice_shelf_dyn(param_file, Time, ISS, CS, G, US, diag, new_
       call initialize_ice_flow_from_file(CS%bed_elev,CS%u_shelf, CS%v_shelf, CS%ground_frac, &
                   G, US, param_file)
       call pass_vector(CS%u_shelf, CS%v_shelf, G%domain, TO_ALL, BGRID_NE)
-      call pass_var(CS%bed_elev, G%domain,CENTER)
+      call pass_var(CS%ground_frac, G%domain)
+      call pass_var(CS%bed_elev, G%domain)
       call update_velocity_masks(CS, G, ISS%hmask, CS%umask, CS%vmask, CS%u_face_mask, CS%v_face_mask)
     endif
   ! Register diagnostics.
@@ -590,8 +601,9 @@ subroutine initialize_ice_shelf_dyn(param_file, Time, ISS, CS, G, US, diag, new_
        'intermediate ocean column thickness passed to ice model', 'm', conversion=US%Z_to_m)
   endif
   call MOM_mesg("MOM_ice_shelf.F90, initialize_ice_shelf: initialize ice velocity.")
-  call update_OD_ffrac_uncoupled(CS, G, ISS%h_shelf(:,:))
-  call ice_shelf_solve_outer(CS, ISS, G, US, CS%u_shelf, CS%v_shelf,CS%taudx_shelf,CS%taudy_shelf, iters, Time)
+  if (new_sim) then
+    call update_OD_ffrac_uncoupled(CS, G, ISS%h_shelf(:,:))
+  endif
 
 end subroutine initialize_ice_shelf_dyn
 
@@ -955,7 +967,7 @@ subroutine ice_shelf_solve_outer(CS, ISS, G, US, u_shlf, v_shlf, taudx, taudy, i
   call CG_action(Au, Av, u_shlf, v_shlf, Phi, Phisub, CS%umask, CS%vmask, ISS%hmask, H_node, &
                  CS%ice_visc, float_cond, CS%bed_elev, CS%basal_traction, &
                  G, US, G%isc-1, G%iec+1, G%jsc-1, G%jec+1, rhoi_rhow)
-  call pass_vector(Au,Av,G%domain,TO_ALL,BGRID_NE)
+  call pass_vector(Au, Av, G%domain, TO_ALL, BGRID_NE)
 
   if (CS%nonlin_solve_err_mode == 1) then
     err_init = 0 ; err_tempu = 0 ; err_tempv = 0
@@ -1012,6 +1024,8 @@ subroutine ice_shelf_solve_outer(CS, ISS, G, US, u_shlf, v_shlf, taudx, taudy, i
                    CS%ice_visc, float_cond, CS%bed_elev, CS%basal_traction, &
                    G, US, G%isc-1, G%iec+1, G%jsc-1, G%jec+1, rhoi_rhow)
 
+    call pass_vector(Au, Av, G%domain, TO_ALL, BGRID_NE)
+
     err_max = 0
 
     if (CS%nonlin_solve_err_mode == 1) then
@@ -1225,7 +1239,7 @@ subroutine ice_shelf_solve_inner(CS, ISS, G, US, u_shlf, v_shlf, taudx, taudy, H
     ! the computational domain - this is their state in the initial iteration
 
 
-    is = isc - cg_halo ; ie = iecq + cg_halo
+    is = iscq - cg_halo ; ie = iecq + cg_halo
     js = jscq - cg_halo ; je = jecq + cg_halo
 
     Au(:,:) = 0 ; Av(:,:) = 0
@@ -2595,10 +2609,9 @@ subroutine calc_shelf_visc(CS, ISS, G, US, u_shlf, v_shlf)
   is = iscq - 1; js = jscq - 1
     i_off = G%idg_offset ; j_off = G%jdg_offset
 
-  allocate(Phi(1:8,1:4,isd:ied,jsd:jed), source=0.0)
+  allocate(Phi(1:8,1:4,isc:iec,jsc:jec), source=0.0)
 
-!  do j=jsc,jec ; do i=isc,iec
-  do j=jscq,jecq ; do i=iscq,iecq
+  do j=jsc,jec ; do i=isc,iec
     call bilinear_shape_fn_grid(G, i, j, Phi(:,:,i,j))
   enddo ; enddo
 

diff --git a/src/ice_shelf/MOM_ice_shelf_initialize.F90 b/src/ice_shelf/MOM_ice_shelf_initialize.F90
@@ -429,31 +429,32 @@ subroutine initialize_ice_flow_from_file(bed_elev,u_shelf, v_shelf,float_cond,&
   filename = trim(inputdir)//trim(vel_file)
   call log_param(PF, mdl, "INPUTDIR/THICKNESS_FILE", filename)
   call get_param(PF, mdl, "ICE_U_VEL_VARNAME", ushelf_varname, &
-                 "The name of the thickness variable in ICE_VELOCITY_FILE.", &
+                 "The name of the u velocity variable in ICE_VELOCITY_FILE.", &
                  default="u_shelf")
   call get_param(PF, mdl, "ICE_V_VEL_VARNAME", vshelf_varname, &
-                 "The name of the thickness variable in ICE_VELOCITY_FILE.", &
+                 "The name of the v velocity variable in ICE_VELOCITY_FILE.", &
                  default="v_shelf")
   call get_param(PF, mdl, "ICE_VISC_VARNAME", ice_visc_varname, &
-                 "The name of the thickness variable in ICE_VELOCITY_FILE.", &
+                 "The name of the ice viscosity variable in ICE_VELOCITY_FILE.", &
                  default="viscosity")
+  call get_param(PF, mdl, "ICE_FLOAT_FRAC_VARNAME", floatfr_varname, &
+                 "The name of the ice float fraction (grounding fraction) variable in ICE_VELOCITY_FILE.", &
+                 default="float_frac")
   call get_param(PF, mdl, "BED_TOPO_FILE", bed_topo_file, &
                  "The file from which the bed elevation is read.", &
                  default="ice_shelf_vel.nc")
   call get_param(PF, mdl, "BED_TOPO_VARNAME", bed_varname, &
-                 "The name of the thickness variable in ICE_INPUT_FILE.", &
+                 "The name of the bed elevation variable in ICE_INPUT_FILE.", &
                  default="depth")
   if (.not.file_exists(filename, G%Domain)) call MOM_error(FATAL, &
        " initialize_ice_shelf_velocity_from_file: Unable to open "//trim(filename))
 
-  floatfr_varname = "float_frac"
-
   call MOM_read_data(filename, trim(ushelf_varname), u_shelf, G%Domain, position=CORNER, scale=US%m_s_to_L_T)
   call MOM_read_data(filename, trim(vshelf_varname), v_shelf, G%Domain, position=CORNER, scale=US%m_s_to_L_T)
   call MOM_read_data(filename, trim(floatfr_varname), float_cond, G%Domain, scale=1.)
 
   filename = trim(inputdir)//trim(bed_topo_file)
-  call MOM_read_data(filename,trim(bed_varname), bed_elev, G%Domain, scale=1.0)
+  call MOM_read_data(filename, trim(bed_varname), bed_elev, G%Domain, scale=US%m_to_Z)
 
 
 end subroutine initialize_ice_flow_from_file