Skip to content

Commit

Permalink
changes to ice_velocity_mask_update, initialization of ice velocity a…
Browse files Browse the repository at this point in the history 
…nd ice thickness for new and restart simulations
  • Loading branch information
@OlgaSergienko
OlgaSergienko committed Aug 25, 2021
1 parent 8ef9596 commit 3ce2efe
Show file tree
Hide file tree
Showing 2 changed files with 34 additions and 34 deletions.
55 changes: 27 additions & 28 deletions src/ice_shelf/MOM_ice_shelf_dynamics.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -277,24 +277,24 @@ subroutine register_ice_shelf_dyn_restarts(G, param_file, CS, restart_CS)
"ice sheet/shelf u-velocity", "m s-1", hor_grid='Bu')
call register_restart_field(CS%v_shelf, "v_shelf", .false., restart_CS, &
"ice sheet/shelf v-velocity", "m s-1", hor_grid='Bu')
call register_restart_field(CS%u_bdry_val, "u_bdry", .false., restart_CS, &
call register_restart_field(CS%u_bdry_val, "u_bdry_val", .false., restart_CS, &
"ice sheet/shelf boundary u-velocity", "m s-1", hor_grid='Bu')
call register_restart_field(CS%v_bdry_val, "v_bdry", .false., restart_CS, &
call register_restart_field(CS%v_bdry_val, "v_bdry_val", .false., restart_CS, &
"ice sheet/shelf boundary v-velocity", "m s-1", hor_grid='Bu')
call register_restart_field(CS%u_face_mask_bdry, "u_bdry_mask", .false., restart_CS, &
call register_restart_field(CS%u_face_mask_bdry, "u_face_mask_bdry", .false., restart_CS, &
"ice sheet/shelf boundary u-mask", "nondim", hor_grid='Bu')
call register_restart_field(CS%v_face_mask_bdry, "v_bdry_mask", .false., restart_CS, &
call register_restart_field(CS%v_face_mask_bdry, "v_face_mask_bdry", .false., restart_CS, &
"ice sheet/shelf boundary v-mask", "nondim", hor_grid='Bu')

call register_restart_field(CS%OD_av, "OD_av", .true., restart_CS, &
"Average open ocean depth in a cell","m")
call register_restart_field(CS%ground_frac, "ground_frac", .true., restart_CS, &
"fractional degree of grounding", "nondim")
call register_restart_field(CS%C_basal_friction, "tau_b_beta", .true., restart_CS, &
call register_restart_field(CS%C_basal_friction, "C_basal_friction", .true., restart_CS, &
"basal sliding coefficients", "Pa (m s-1)^n_sliding")
call register_restart_field(CS%AGlen_visc, "A_Glen", .true., restart_CS, &
call register_restart_field(CS%AGlen_visc, "AGlen_visc", .true., restart_CS, &
"ice-stiffness parameter", "Pa-3 s-1")
call register_restart_field(CS%h_bdry_val, "h_bdry", .false., restart_CS, &
call register_restart_field(CS%h_bdry_val, "h_bdry_val", .false., restart_CS, &
"ice thickness at the boundary","m")
endif

Expand Down Expand Up @@ -503,6 +503,7 @@ subroutine initialize_ice_shelf_dyn(param_file, Time, ISS, CS, G, US, diag, new_
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)
endif

Expand Down Expand Up @@ -533,6 +534,7 @@ subroutine initialize_ice_shelf_dyn(param_file, Time, ISS, CS, G, US, diag, new_
endif

! 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)

Expand All @@ -556,24 +558,23 @@ subroutine initialize_ice_shelf_dyn(param_file, Time, ISS, CS, G, US, diag, new_
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 update_velocity_masks(CS, G, ISS%hmask, CS%umask, CS%vmask, CS%u_face_mask, CS%v_face_mask)

endif
! Register diagnostics.
CS%id_u_shelf = register_diag_field('ice_shelf_model','u_shelf',CS%diag%axesB1, Time, &
'x-velocity of ice', 'm yr-1', conversion=365.0*86400.0*US%L_T_to_m_s)
CS%id_v_shelf = register_diag_field('ice_shelf_model','v_shelf',CS%diag%axesB1, Time, &
'y-velocity of ice', 'm yr-1', conversion=365.0*86400.0*US%L_T_to_m_s)
! I think that the conversion factors for the next two diagnostics are wrong. - RWH
CS%id_taudx_shelf = register_diag_field('ice_shelf_model','taudx_shelf',CS%diag%axesB1, Time, &
'x-driving stress of ice', 'kPa', conversion=1.e-9*US%L_T_to_m_s)
'x-driving stress of ice', 'kPa', conversion=1.e-9*US%RL2_T2_to_Pa)
CS%id_taudy_shelf = register_diag_field('ice_shelf_model','taudy_shelf',CS%diag%axesB1, Time, &
'y-driving stress of ice', 'kPa', conversion=1.e-9*US%L_T_to_m_s)
'y-driving stress of ice', 'kPa', conversion=1.e-9*US%RL2_T2_to_Pa)
CS%id_u_mask = register_diag_field('ice_shelf_model','u_mask',CS%diag%axesB1, Time, &
'mask for u-nodes', 'none')
CS%id_v_mask = register_diag_field('ice_shelf_model','v_mask',CS%diag%axesB1, Time, &
'mask for v-nodes', 'none')
CS%id_ground_frac = register_diag_field('ice_shelf_model','ice_ground_frac',CS%diag%axesT1, Time, &
'fraction of cell that is grounded', 'none')
!
CS%id_col_thick = register_diag_field('ice_shelf_model','col_thick',CS%diag%axesT1, Time, &
'ocean column thickness passed to ice model', 'm', conversion=US%Z_to_m)
CS%id_visc_shelf = register_diag_field('ice_shelf_model','ice_visc',CS%diag%axesT1, Time, &
Expand All @@ -582,12 +583,10 @@ subroutine initialize_ice_shelf_dyn(param_file, Time, ISS, CS, G, US, diag, new_
'taub', 'Pa yr m-1', conversion=1e-6*US%Z_to_m)
CS%id_OD_av = register_diag_field('ice_shelf_model','OD_av',CS%diag%axesT1, Time, &
'intermediate ocean column thickness passed to ice model', 'm', conversion=US%Z_to_m)
if (new_sim) then
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)
endif
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)

end subroutine initialize_ice_shelf_dyn

Expand Down Expand Up @@ -960,7 +959,7 @@ subroutine ice_shelf_solve_outer(CS, ISS, G, US, u_shlf, v_shlf,taudx,taudy, ite

!! begin loop

do iter=1,400
do iter=1,100

call ice_shelf_solve_inner(CS, ISS, G, US, u_shlf, v_shlf, taudx, taudy, H_node, float_cond, &
ISS%hmask, conv_flag, iters, time, Phi, Phisub)
Expand Down Expand Up @@ -1043,16 +1042,15 @@ subroutine ice_shelf_solve_outer(CS, ISS, G, US, u_shlf, v_shlf,taudx,taudy, ite
call MOM_mesg(mesg, 5)

if (err_max <= CS%nonlinear_tolerance * err_init) then
write(mesg,*) "ice_shelf_solve_outer: nonlinear fractional residual = ", err_max/err_init
call MOM_mesg(mesg)
write(mesg,*) "ice_shelf_solve_outer: exiting nonlinear solve after ",iter," iterations"
! call MOM_mesg(mesg, 5)
call MOM_mesg(mesg)
exit
endif

enddo

write(mesg,*) "ice_shelf_solve_outer: nonlinear fractional residual = ", err_max/err_init
call MOM_mesg(mesg)
write(mesg,*) "ice_shelf_solve_outer: exiting nonlinear solve after ",iter," iterations"
call MOM_mesg(mesg)
deallocate(Phi)
deallocate(Phisub)

Expand Down Expand Up @@ -1086,6 +1084,7 @@ subroutine ice_shelf_solve_inner(CS, ISS, G, US, u_shlf, v_shlf, taudx, taudy, H
!! iterations have converged to the specified tolerance
integer, intent(out) :: iters !< The number of iterations used in the solver.
type(time_type), intent(in) :: Time !< The current model time
character(len=160) :: mesg ! The text of an error message
real, dimension(8,4,SZDI_(G),SZDJ_(G)), &
intent(in) :: Phi !< The gradients of bilinear basis elements at Gaussian
!! quadrature points surrounding the cell vertices [L-1 ~> m-1].
Expand Down Expand Up @@ -2586,7 +2585,7 @@ subroutine calc_shelf_visc(CS, ISS, G, US, u_shlf, v_shlf)
(u_shlf(I,J-1) + (u_shlf(I-1,J-1) + u_shlf(I+1,J-1)))) / (3*G%dyT(i,j))
vy = ((v_shlf(I,J) + (v_shlf(I-1,J)+ v_shlf(I+1,J))) - &
(v_shlf(I,J-1) + (v_shlf(I-1,J-1)+ v_shlf(I+1,J-1)))) / (3*G%dyT(i,j))
! CS%ice_visc(i,j) =1e14*(G%areaT(i,j) * ISS%h_shelf(i,j)) ! constant viscocity for debugging
! CS%ice_visc(i,j) =1e15*(G%areaT(i,j) * ISS%h_shelf(i,j)) ! constant viscocity for debugging
CS%ice_visc(i,j) = 0.5 * Visc_coef * (G%areaT(i,j) * ISS%h_shelf(i,j)) * &
(US%s_to_T**2 * (ux**2 + vy**2 + ux*vy + 0.25*(uy+vx)**2 + eps_min**2))**((1.-n_g)/(2.*n_g))
endif
Expand Down Expand Up @@ -2938,7 +2937,7 @@ subroutine update_velocity_masks(CS, G, hmask, umask, vmask, u_face_mask, v_face
do j=js,G%jed
do i=is,G%ied

if (hmask(i,j) == 1) then
if ((hmask(i,j) == 1) .OR. (hmask(i,j) == 3)) then

umask(I,j) = 1.
vmask(I,j) = 1.
Expand All @@ -2947,10 +2946,10 @@ subroutine update_velocity_masks(CS, G, hmask, umask, vmask, u_face_mask, v_face

select case (int(CS%u_face_mask_bdry(I-1+k,j)))
case (3)
! vmask(I-1+k,J-1)=0.
vmask(I-1+k,J-1)=3.
u_face_mask(I-1+k,j)=3.
umask(I-1+k,J)=3.
!vmask(I-1+k,J)=0.
vmask(I-1+k,J)=3.
vmask(I-1+k,J)=3.
case (2)
u_face_mask(I-1+k,j)=2.
Expand All @@ -2973,9 +2972,9 @@ subroutine update_velocity_masks(CS, G, hmask, umask, vmask, u_face_mask, v_face
select case (int(CS%v_face_mask_bdry(i,J-1+k)))
case (3)
vmask(I-1,J-1+k)=3.
umask(I-1,J-1+k)=0.
umask(I-1,J-1+k)=3.
vmask(I,J-1+k)=3.
umask(I,J-1+k)=0.
umask(I,J-1+k)=3.
v_face_mask(i,J-1+k)=3.
case (2)
v_face_mask(i,J-1+k)=2.
Expand Down
13 changes: 7 additions & 6 deletions src/ice_shelf/MOM_ice_shelf_initialize.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -101,7 +101,7 @@ subroutine initialize_ice_thickness_from_file(h_shelf, area_shelf_h, hmask, G, U
! 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,thickness_file,inputdir ! Strings for file/path
character(len=200) :: thickness_varname, area_varname ! Variable name in file
character(len=200) :: thickness_varname, area_varname, hmask_varname ! Variable name in file
character(len=40) :: mdl = "initialize_ice_thickness_from_file" ! This subroutine's name.
integer :: i, j, isc, jsc, iec, jec
real :: len_sidestress, mask, udh
Expand All @@ -125,22 +125,22 @@ subroutine initialize_ice_thickness_from_file(h_shelf, area_shelf_h, hmask, G, U
call get_param(PF, mdl, "ICE_AREA_VARNAME", area_varname, &
"The name of the area variable in ICE_THICKNESS_FILE.", &
default="area_shelf_h")

hmask_varname="h_mask"
if (.not.file_exists(filename, G%Domain)) call MOM_error(FATAL, &
" initialize_topography_from_file: Unable to open "//trim(filename))
call MOM_read_data(filename, trim(thickness_varname), h_shelf, G%Domain, scale=US%m_to_Z)
call MOM_read_data(filename,trim(area_varname), area_shelf_h, G%Domain, scale=US%m_to_L**2)

call MOM_read_data(filename, trim(hmask_varname), hmask, G%Domain)
isc = G%isc ; jsc = G%jsc ; iec = G%iec ; jec = G%jec

if (len_sidestress > 0.) then

This comment has been minimized.

Sign in to view

Copy link
@MJHarrison-GFDL

MJHarrison-GFDL Dec 1, 2021

Contributor

This conditioinal block changes answers in certain ice shelf configurations where DYNAMIC_SHELF_MASS = False and LEN_SIDE_STRESS = 0.0 . Was there a reason to apply this where DYNAMIC_SHELF_MASS = True ? Also, is hmask now required? This changes answers where previously, hmask was derived from area_shelf.

This comment has been minimized.

Sign in to view

Copy link
@OlgaSergienko

OlgaSergienko Dec 1, 2021

Author Contributor

This is for idealized configurations, where tapering was applied. For realistic configurations len_sidestress =0
do j=jsc,jec
do i=isc,iec

! taper ice shelf in area where there is no sidestress -
! but do not interfere with hmask

if ((G%geoLonCv(i,j) > len_sidestress).and. &
(len_sidestress > 0.)) then
if (G%geoLonCv(i,j) > len_sidestress) then
udh = exp(-(G%geoLonCv(i,j)-len_sidestress)/5.0) * h_shelf(i,j)
if (udh <= 25.0) then
h_shelf(i,j) = 0.0
Expand All @@ -154,6 +154,7 @@ subroutine initialize_ice_thickness_from_file(h_shelf, area_shelf_h, hmask, G, U

if (area_shelf_h (i,j) >= G%areaT(i,j)) then
hmask(i,j) = 1.
area_shelf_h(i,j)=G%areaT(i,j)
elseif (area_shelf_h (i,j) == 0.0) then
hmask(i,j) = 0.
elseif ((area_shelf_h(i,j) > 0) .and. (area_shelf_h(i,j) <= G%areaT(i,j))) then
Expand All @@ -163,7 +164,7 @@ subroutine initialize_ice_thickness_from_file(h_shelf, area_shelf_h, hmask, G, U
endif
enddo
enddo

endif
end subroutine initialize_ice_thickness_from_file

!> Initialize ice shelf thickness for a channel configuration
Expand Down

0 comments on commit 3ce2efe

Please sign in to comment.