Skip to content

Commit

Permalink
A new subrutine to initialize the bed elevation under the ice sheet f…
Browse files Browse the repository at this point in the history 
…rom an input file separetely from other fields; fix to a bug Claire has found.
  • Loading branch information
@OlgaSergienko
OlgaSergienko committed Sep 6, 2023
1 parent e9d1cb5 commit f14ec1e
Show file tree
Hide file tree
Showing 3 changed files with 89 additions and 32 deletions.
2 changes: 1 addition & 1 deletion src/ice_shelf/MOM_ice_shelf.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -637,7 +637,7 @@ subroutine shelf_calc_flux(sfc_state_in, fluxes_in, Time, time_step_in, CS)
Sbdry(i,j) = Sbdry_it
endif ! Sb_min_set

Sbdry(i,j) = Sbdry_it
! Sbdry(i,j) = Sbdry_it
endif ! CS%find_salt_root

enddo !it1
Expand Down
61 changes: 43 additions & 18 deletions src/ice_shelf/MOM_ice_shelf_dynamics.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -25,7 +25,7 @@ module MOM_ice_shelf_dynamics
use MOM_checksums, only : hchksum, qchksum
use MOM_ice_shelf_initialize, only : initialize_ice_shelf_boundary_channel,initialize_ice_flow_from_file
use MOM_ice_shelf_initialize, only : initialize_ice_shelf_boundary_from_file,initialize_ice_C_basal_friction
use MOM_ice_shelf_initialize, only : initialize_ice_AGlen
use MOM_ice_shelf_initialize, only : initialize_ice_AGlen, initialize_ice_bed_elevation
implicit none ; private

#include <MOM_memory.h>
Expand Down Expand Up @@ -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_elev", .false., restart_CS, &
"bed elevation", "m", conversion=US%Z_to_m)
endif

end subroutine register_ice_shelf_dyn_restarts
Expand Down Expand Up @@ -481,7 +483,7 @@ subroutine initialize_ice_shelf_dyn(param_file, Time, ISS, CS, G, US, diag, new_

! Take additional initialization steps, for example of dependent variables.
if (active_shelf_dynamics .and. .not.new_sim) then

call initialize_ice_bed_elevation(CS%bed_elev, G, US, param_file)
! this is unfortunately necessary; if grid is not symmetric the boundary values
! of u and v are otherwise not set till the end of the first linear solve, and so
! viscosity is not calculated correctly.
Expand All @@ -504,12 +506,31 @@ subroutine initialize_ice_shelf_dyn(param_file, Time, ISS, CS, G, US, diag, new_
enddo ; enddo
endif

call pass_var(CS%OD_av,G%domain)
call pass_var(CS%ground_frac,G%domain)
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_vector(CS%u_shelf, CS%v_shelf, G%domain, TO_ALL, BGRID_NE)

call pass_var(CS%OD_av,G%domain, complete=.false.)
call pass_var(CS%ground_frac,G%domain, complete=.false.)
call pass_var(CS%ice_visc,G%domain, complete=.false.)
call pass_var(CS%basal_traction, G%domain, complete=.false.)
call pass_var(CS%AGlen_visc, G%domain, complete=.false.)
call pass_var(CS%bed_elev, G%domain, complete=.false.)
call pass_var(CS%C_basal_friction, G%domain, complete=.false.)
call pass_var(CS%h_bdry_val, G%domain, complete=.false.)
call pass_var(CS%thickness_bdry_val, G%domain, complete=.true.)

call pass_vector(CS%u_shelf, CS%v_shelf, G%domain, TO_ALL, BGRID_NE, complete=.false.)
call pass_vector(CS%u_bdry_val, CS%v_bdry_val, G%domain, TO_ALL, BGRID_NE, complete=.false.)
call pass_vector(CS%u_face_mask_bdry, CS%v_face_mask_bdry, G%domain, TO_ALL, BGRID_NE, complete=.true.)
call update_velocity_masks(CS, G, ISS%hmask, CS%umask, CS%vmask, CS%u_face_mask, CS%v_face_mask)
! call pass_var(CS%OD_av,G%domain)
! call pass_var(CS%ground_frac,G%domain)
! 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_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, complete=.false.)
! call pass_vector(CS%u_face_mask_bdry, CS%v_face_mask_bdry, G%domain, TO_ALL, BGRID_NE, complete=.true.)
! 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
Expand Down Expand Up @@ -541,27 +562,31 @@ subroutine initialize_ice_shelf_dyn(param_file, Time, ISS, CS, G, US, diag, new_
! initialize basal friction coefficients
if (new_sim) then
call initialize_ice_C_basal_friction(CS%C_basal_friction, G, US, param_file)
call pass_var(CS%C_basal_friction, G%domain)
call pass_var(CS%C_basal_friction, G%domain, complete=.false.)

! initialize ice-stiffness AGlen
call initialize_ice_AGlen(CS%AGlen_visc, G, US, param_file)
call pass_var(CS%AGlen_visc, G%domain)
call pass_var(CS%AGlen_visc, G%domain, complete=.false.)

!initialize boundary conditions
call initialize_ice_shelf_boundary_from_file(CS%u_face_mask_bdry, CS%v_face_mask_bdry, &
CS%u_bdry_val, CS%v_bdry_val, CS%umask, CS%vmask, CS%h_bdry_val, &
CS%thickness_bdry_val, ISS%hmask, ISS%h_shelf, G, US, param_file )
call pass_var(ISS%hmask, 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_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 pass_var(ISS%hmask, G%domain, complete=.false.)
call pass_var(CS%h_bdry_val, G%domain, complete=.false.)
call pass_var(CS%thickness_bdry_val, G%domain, complete=.true.)
call pass_vector(CS%u_bdry_val, CS%v_bdry_val, G%domain, TO_ALL, BGRID_NE, complete=.false.)
call pass_vector(CS%u_face_mask_bdry, CS%v_face_mask_bdry, G%domain, TO_ALL, BGRID_NE, complete=.false.)

!initialize ice flow characteristic (velocities, bed elevation under the grounded part, etc) from file
call initialize_ice_flow_from_file(CS%bed_elev,CS%u_shelf, CS%v_shelf, CS%ground_frac, &
! call initialize_ice_flow_from_file(CS%bed_elev,CS%u_shelf, CS%v_shelf, CS%ground_frac, &
! G, US, param_file)
call initialize_ice_flow_from_file(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_vector(CS%u_shelf, CS%v_shelf, G%domain, TO_ALL, BGRID_NE, complete=.true.)
call pass_var(CS%ground_frac,G%domain, complete=.false.)
call initialize_ice_bed_elevation(CS%bed_elev, G, US, param_file )
call pass_var(CS%bed_elev, G%domain, complete=.true.)
call update_velocity_masks(CS, G, ISS%hmask, CS%umask, CS%vmask, CS%u_face_mask, CS%v_face_mask)
endif
! Register diagnostics.
Expand Down
58 changes: 45 additions & 13 deletions src/ice_shelf/MOM_ice_shelf_initialize.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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_bed_elevation
! 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
Expand Down Expand Up @@ -390,13 +391,15 @@ end subroutine initialize_ice_shelf_boundary_channel


!> Initialize ice shelf flow from file
subroutine initialize_ice_flow_from_file(bed_elev,u_shelf, v_shelf,float_cond,&
subroutine initialize_ice_flow_from_file(u_shelf, v_shelf,float_cond,&
G, US, PF)
!subroutine initialize_ice_flow_from_file(bed_elev,u_shelf, v_shelf,float_cond,&
! G, US, PF)
!subroutine initialize_ice_flow_from_file(bed_elev,u_shelf, v_shelf,float_cond,ice_visc,&
! G, US, PF)
type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure
real, dimension(SZDI_(G),SZDJ_(G)), &
intent(inout) :: bed_elev !< The bed elevation [Z ~> m].
! real, dimension(SZDI_(G),SZDJ_(G)), &
! intent(inout) :: bed_elev !< The bed elevation [Z ~> m].
real, dimension(SZIB_(G),SZJB_(G)), &
intent(inout) :: u_shelf !< The zonal ice shelf velocity [L T-1 ~> m s-1].
real, dimension(SZIB_(G),SZJB_(G)), &
Expand All @@ -409,9 +412,9 @@ subroutine initialize_ice_flow_from_file(bed_elev,u_shelf, v_shelf,float_cond,&

! This subroutine reads ice thickness and area from a file and puts it into
! h_shelf [Z ~> m] and area_shelf_h [L2 ~> m2] (and dimensionless) and updates hmask
character(len=200) :: filename,vel_file,inputdir,bed_topo_file ! Strings for file/path
character(len=200) :: filename,vel_file,inputdir,i!bed_topo_file ! Strings for file/path
character(len=200) :: ushelf_varname, vshelf_varname, &
ice_visc_varname, floatfr_varname, bed_varname ! Variable name in file
ice_visc_varname, floatfr_varname!, bed_varname ! Variable name in file
character(len=40) :: mdl = "initialize_ice_velocity_from_file" ! This subroutine's name.
real :: len_sidestress

Expand All @@ -437,12 +440,12 @@ subroutine initialize_ice_flow_from_file(bed_elev,u_shelf, v_shelf,float_cond,&
call get_param(PF, mdl, "ICE_VISC_VARNAME", ice_visc_varname, &
"The name of the thickness variable in ICE_VELOCITY_FILE.", &
default="viscosity")
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.", &
default="depth")
! 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.", &
! default="depth")
if (.not.file_exists(filename, G%Domain)) call MOM_error(FATAL, &
" initialize_ice_shelf_velocity_from_file: Unable to open "//trim(filename))

Expand All @@ -452,8 +455,8 @@ subroutine initialize_ice_flow_from_file(bed_elev,u_shelf, v_shelf,float_cond,&
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)
! filename = trim(inputdir)//trim(bed_topo_file)
! call MOM_read_data(filename,trim(bed_varname), bed_elev, G%Domain, scale=1.0)


end subroutine initialize_ice_flow_from_file
Expand Down Expand Up @@ -657,4 +660,33 @@ subroutine initialize_ice_AGlen(AGlen, G, US, PF)

endif
end subroutine

!> Initialize bed elevation B
subroutine initialize_ice_bed_elevation(bed_elev, G, US, PF)
type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure
real, dimension(SZDI_(G),SZDJ_(G)), &
intent(inout) :: bed_elev !< The bed elevation under the ice sheet [m]
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 :: A_Glen ! Ice-stiffness parameter, often in [Pa-3 s-1]
character(len=40) :: mdl = "initialize_ice_bed_elev" ! This subroutine's name.
character(len=200) :: config
character(len=200) :: bed_varname
character(len=200) :: inputdir, filename, bed_topo_file

call get_param(PF, mdl, "INPUTDIR", inputdir, default=".")
inputdir = slasher(inputdir)
call get_param(PF, mdl, "BED_TOPO_FILE", bed_topo_file, &
"The file from which the bed elevation is read.", &
default="IS_input.nc")
call get_param(PF, mdl, "BED_TOPO_VARNAME", bed_varname, &
"The name of the bed elevation variable in BED_TOPO_FILE.", &
default="bed_elev")
filename = trim(inputdir)//trim(bed_topo_file)
if (.not.file_exists(filename, G%Domain)) call MOM_error(FATAL, &
" initialize_ice_bed_elevation_from_file: Unable to open "//trim(filename))
call MOM_read_data(filename,trim(bed_varname), bed_elev, G%Domain, scale=1.0)

end subroutine
end module MOM_ice_shelf_initialize

0 comments on commit f14ec1e

Please sign in to comment.