Skip to content

Commit

Permalink
Merge pull request #1479 from OlgaSergienko/ice_dynamics
Browse files Browse the repository at this point in the history 
Ice dynamics
  • Loading branch information
@marshallward
marshallward authored Sep 1, 2021
2 parents f98b76d + d234258 commit a7d8e3a
Show file tree
Hide file tree
Showing 2 changed files with 104 additions and 104 deletions.
159 changes: 79 additions & 80 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,47 +534,47 @@ subroutine initialize_ice_shelf_dyn(param_file, Time, ISS, CS, G, US, diag, new_
endif

! initialize basal friction coefficients
call initialize_ice_C_basal_friction(CS%C_basal_friction, G, US, param_file)
call pass_var(CS%C_basal_friction, G%domain)
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)

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

!initialize boundary conditions
call initialize_ice_shelf_boundary_from_file(CS%u_face_mask_bdry, CS%v_face_mask_bdry, &
!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)

!initialize ice flow velocities from file
call initialize_ice_flow_from_file(CS%bed_elev,CS%u_shelf, CS%v_shelf,CS%ground_frac, ISS%hmask,ISS%h_shelf, &
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)

!initialize ice flow velocities from file
call initialize_ice_flow_from_file(CS%bed_elev,CS%u_shelf, CS%v_shelf,CS%ground_frac, ISS%hmask,ISS%h_shelf, &
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 update_velocity_masks(CS, G, ISS%hmask, CS%umask, CS%vmask, CS%u_face_mask, CS%v_face_mask)

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 @@ -683,46 +682,46 @@ subroutine update_ice_shelf(CS, ISS, G, US, time_step, Time, ocean_mass, coupled
coupled_GL = .false.
if (present(ocean_mass) .and. present(coupled_grounding)) coupled_GL = coupled_grounding
!
call ice_shelf_advect(CS, ISS, G, time_step, Time)
CS%elapsed_velocity_time = CS%elapsed_velocity_time + time_step
if (CS%elapsed_velocity_time >= CS%velocity_update_time_step) update_ice_vel = .true.

if (coupled_GL) then
call update_OD_ffrac(CS, G, US, ocean_mass, update_ice_vel)
elseif (update_ice_vel) then
call update_OD_ffrac_uncoupled(CS, G, ISS%h_shelf(:,:))
endif
call ice_shelf_advect(CS, ISS, G, time_step, Time)
CS%elapsed_velocity_time = CS%elapsed_velocity_time + time_step
if (CS%elapsed_velocity_time >= CS%velocity_update_time_step) update_ice_vel = .true.

if (coupled_GL) then
call update_OD_ffrac(CS, G, US, ocean_mass, update_ice_vel)
elseif (update_ice_vel) then
call update_OD_ffrac_uncoupled(CS, G, ISS%h_shelf(:,:))
endif

if (update_ice_vel) then
call ice_shelf_solve_outer(CS, ISS, G, US, CS%u_shelf, CS%v_shelf,CS%taudx_shelf,CS%taudy_shelf, iters, Time)
endif

if (update_ice_vel) then
call ice_shelf_solve_outer(CS, ISS, G, US, CS%u_shelf, CS%v_shelf,CS%taudx_shelf,CS%taudy_shelf, iters, Time)
endif

! call ice_shelf_temp(CS, ISS, G, US, time_step, ISS%water_flux, Time)

if (update_ice_vel) 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)
if (CS%id_v_shelf > 0) call post_data(CS%id_v_shelf, CS%v_shelf, CS%diag)
if (update_ice_vel) 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)
if (CS%id_v_shelf > 0) call post_data(CS%id_v_shelf, CS%v_shelf, CS%diag)
! if (CS%id_t_shelf > 0) call post_data(CS%id_t_shelf,CS%t_shelf,CS%diag)
if (CS%id_taudx_shelf > 0) call post_data(CS%id_taudx_shelf, CS%taudx_shelf, CS%diag)
if (CS%id_taudy_shelf > 0) call post_data(CS%id_taudy_shelf, CS%taudy_shelf, CS%diag)
if (CS%id_ground_frac > 0) call post_data(CS%id_ground_frac, CS%ground_frac,CS%diag)
if (CS%id_OD_av >0) call post_data(CS%id_OD_av, CS%OD_av,CS%diag)
if (CS%id_visc_shelf > 0) call post_data(CS%id_visc_shelf, CS%ice_visc,CS%diag)
if (CS%id_taub > 0) call post_data(CS%id_taub, CS%basal_traction,CS%diag)
if (CS%id_taudx_shelf > 0) call post_data(CS%id_taudx_shelf, CS%taudx_shelf, CS%diag)
if (CS%id_taudy_shelf > 0) call post_data(CS%id_taudy_shelf, CS%taudy_shelf, CS%diag)
if (CS%id_ground_frac > 0) call post_data(CS%id_ground_frac, CS%ground_frac,CS%diag)
if (CS%id_OD_av >0) call post_data(CS%id_OD_av, CS%OD_av,CS%diag)
if (CS%id_visc_shelf > 0) call post_data(CS%id_visc_shelf, CS%ice_visc,CS%diag)
if (CS%id_taub > 0) call post_data(CS%id_taub, CS%basal_traction,CS%diag)
!!
if (CS%id_u_mask > 0) call post_data(CS%id_u_mask,CS%umask,CS%diag)
if (CS%id_v_mask > 0) call post_data(CS%id_v_mask,CS%vmask,CS%diag)
if (CS%id_ufb_mask > 0) call post_data(CS%id_ufb_mask,CS%u_face_mask_bdry,CS%diag)
if (CS%id_vfb_mask > 0) call post_data(CS%id_vfb_mask,CS%v_face_mask_bdry,CS%diag)
if (CS%id_u_mask > 0) call post_data(CS%id_u_mask,CS%umask,CS%diag)
if (CS%id_v_mask > 0) call post_data(CS%id_v_mask,CS%vmask,CS%diag)
if (CS%id_ufb_mask > 0) call post_data(CS%id_ufb_mask,CS%u_face_mask_bdry,CS%diag)
if (CS%id_vfb_mask > 0) call post_data(CS%id_vfb_mask,CS%v_face_mask_bdry,CS%diag)
! if (CS%id_t_mask > 0) call post_data(CS%id_t_mask,CS%tmask,CS%diag)

call disable_averaging(CS%diag)
call disable_averaging(CS%diag)

CS%elapsed_velocity_time = 0.0
endif
CS%elapsed_velocity_time = 0.0
endif

end subroutine update_ice_shelf

Expand Down Expand Up @@ -869,13 +868,13 @@ subroutine ice_shelf_solve_outer(CS, ISS, G, US, u_shlf, v_shlf,taudx,taudy, ite
CS%ground_frac(:,:) = 0.0
allocate(Phisub(nsub,nsub,2,2,2,2)) ; Phisub(:,:,:,:,:,:) = 0.0

do j=G%jsc,G%jec
do i=G%isc,G%iec
if (rhoi_rhow * ISS%h_shelf(i,j) - G%bathyT(i,j) > 0) then
do j=G%jsc,G%jec
do i=G%isc,G%iec
if (rhoi_rhow * ISS%h_shelf(i,j) - G%bathyT(i,j) > 0) then
float_cond(i,j) = 1.0
CS%ground_frac(i,j) = 1.0
endif
enddo
endif
enddo
enddo

call calc_shelf_driving_stress(CS, ISS, G, US, taudx, taudy, CS%OD_av)
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
49 changes: 25 additions & 24 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,45 +125,46 @@ 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

do j=jsc,jec
do i=isc,iec
if (len_sidestress > 0.) then
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
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
area_shelf_h(i,j) = 0.0
else
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
area_shelf_h(i,j) = 0.0
else
h_shelf(i,j) = udh
endif
endif
endif

! update thickness mask

if (area_shelf_h (i,j) >= G%areaT(i,j)) then
hmask(i,j) = 1.
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
hmask(i,j) = 2.
else
call MOM_error(FATAL,mdl// " AREA IN CELL OUT OF RANGE")
endif
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
hmask(i,j) = 2.
else
call MOM_error(FATAL,mdl// " AREA IN CELL OUT OF RANGE")
endif
enddo
enddo
enddo

endif
end subroutine initialize_ice_thickness_from_file

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

0 comments on commit a7d8e3a

Please sign in to comment.