Skip to content

Commit

Permalink
Ice dynamics (#228)
Browse files Browse the repository at this point in the history
In this PR an option is added to use ice viscosity computed from the observed surface velocity, computed by the model and use a constant value (for debugging purposes). A new (char) parameter "ICE_VISCOSITY_COMPUTE" is introduced; its values can be "MODEL" (the ice viscosity computed by the model); "OBS" the ice viscosity is computed at the preprocessing step and read from a file (its name is defined by the parameter "ICE_STIFFNESS_FILE") into a variable with a name defined by "A_GLEN_VARNAME" parameter; "CONSANT" is a constant value defined by a parameter "A_GLEN". These changes are in MOM_ice_shelf_dynamics.F90. Minor changes are done to MOM_ice_shelf_initialize.F90 to correct units, scales.
  • Loading branch information
OlgaSergienko authored Nov 11, 2022
1 parent 9d5a320 commit f92a4ab
Show file tree
Hide file tree
Showing 2 changed files with 45 additions and 30 deletions.
51 changes: 33 additions & 18 deletions src/ice_shelf/MOM_ice_shelf_dynamics.F90
Original file line number Diff line number Diff line change
Expand Up @@ -77,7 +77,8 @@ module MOM_ice_shelf_dynamics
real, pointer, dimension(:,:) :: t_shelf => NULL() !< Vertically integrated temperature in the ice shelf/stream,
!! on corner-points (B grid) [C ~> degC]
real, pointer, dimension(:,:) :: tmask => NULL() !< A mask on tracer points that is 1 where there is ice.
real, pointer, dimension(:,:) :: ice_visc => NULL() !< Glen's law ice viscosity, often in [R L4 Z T-1 ~> kg m2 s-1].
real, pointer, dimension(:,:) :: ice_visc => NULL() !< Glen's law ice viscosity (Pa s),
!! in [R L2 T-1 ~> kg m-1 s-1].
real, pointer, dimension(:,:) :: AGlen_visc => NULL() !< Ice-stiffness parameter in Glen's law ice viscosity,
!! often in [kg-1/3 m-1/3 s-1].
real, pointer, dimension(:,:) :: thickness_bdry_val => NULL() !< The ice thickness at an inflowing boundary [Z ~> m].
Expand All @@ -94,7 +95,7 @@ module MOM_ice_shelf_dynamics
!! Sign convention: positive below sea-level, negative above.

real, pointer, dimension(:,:) :: basal_traction => NULL() !< The area integrated nonlinear part of "linearized"
!! basal stress [R Z L2 T-1 ~> kg s-1].
!! basal stress (Pa) [R L2 T-2 ~> Pa].
!! The exact form depends on basal law exponent and/or whether flow is "hybridized" a la Goldberg 2011
real, pointer, dimension(:,:) :: C_basal_friction => NULL()!< Coefficient in sliding law tau_b = C u^(n_basal_fric),
!! units= Pa (m yr-1)-(n_basal_fric)
Expand All @@ -117,6 +118,9 @@ module MOM_ice_shelf_dynamics
real :: g_Earth !< The gravitational acceleration [L2 Z-1 T-2 ~> m s-2].
real :: density_ice !< A typical density of ice [R ~> kg m-3].

character(len=40) :: ice_viscosity_compute !< Specifies whether the ice viscosity is computed internally
!! according to Glen's flow law; is constant (for debugging purposes)
!! or using observed strain rates and read from a file
logical :: GL_regularize !< Specifies whether to regularize the floatation condition
!! at the grounding line as in Goldberg Holland Schoof 2009
integer :: n_sub_regularize
Expand Down Expand Up @@ -261,7 +265,7 @@ subroutine register_ice_shelf_dyn_restarts(G, US, param_file, CS, restart_CS)
allocate( CS%t_shelf(isd:ied,jsd:jed), source=-10.0*US%degC_to_C ) ! [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 )
allocate( CS%basal_traction(isd:ied,jsd:jed), source=0.0 ) ! [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 )
Expand Down Expand Up @@ -429,6 +433,10 @@ subroutine initialize_ice_shelf_dyn(param_file, Time, ISS, CS, G, US, diag, new_
call get_param(param_file, mdl, "CALVE_TO_MASK", CS%calve_to_mask, &
"If true, do not allow an ice shelf where prohibited by a mask.", &
default=.false.)
call get_param(param_file, mdl, "ICE_VISCOSITY_COMPUTE", CS%ice_viscosity_compute, &
"If MODEL, compute ice viscosity internally, if OBS read from a file,"//&
"if CONSTANT a constant value (for debugging).", &
default="MODEL")

endif
call get_param(param_file, mdl, "MIN_THICKNESS_SIMPLE_CALVE", CS%min_thickness_simple_calve, &
Expand Down Expand Up @@ -581,7 +589,7 @@ subroutine initialize_ice_shelf_dyn(param_file, Time, ISS, CS, G, US, diag, new_
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, &
'vi-viscosity', 'Pa s-1 m', conversion=US%RL2_T2_to_Pa*US%L_T_to_m_s) !vertically integrated viscosity
'vi-viscosity', 'Pa m s', conversion=US%RL2_T2_to_Pa*US%Z_to_m*US%T_to_s) !vertically integrated viscosity
CS%id_taub = register_diag_field('ice_shelf_model','taub_beta',CS%diag%axesT1, Time, &
'taub', 'MPa', conversion=1e-6*US%RL2_T2_to_Pa)
CS%id_OD_av = register_diag_field('ice_shelf_model','OD_av',CS%diag%axesT1, Time, &
Expand Down Expand Up @@ -1961,13 +1969,12 @@ subroutine calc_shelf_driving_stress(CS, ISS, G, US, taudx, taudy, OD)
taudy(I,J) = taudy(I,J) - .25 * rho * grav * ISS%h_shelf(i,j) * sy * G%areaT(i,j)
endif
if (CS%ground_frac(i,j) == 1) then
! neumann_val = (.5 * grav * (rho * ISS%h_shelf(i,j)**2 - rhow * CS%bed_elev(i,j)**2))
neumann_val = .5 * grav * (1-rho/rhow) * rho * ISS%h_shelf(i,j)**2
neumann_val = (.5 * grav * (rho * ISS%h_shelf(i,j)**2 - rhow * CS%bed_elev(i,j)**2))
else
neumann_val = (.5 * grav * (1-rho/rhow) * rho * ISS%h_shelf(i,j)**2)
endif

if ((CS%u_face_mask(I-1,j) == 2) .OR. (ISS%hmask(i-1,j) == 0) .OR. (ISS%hmask(i-1,j) == 2) ) then
if ((CS%u_face_mask_bdry(I-1,j) == 2) .OR. (ISS%hmask(i-1,j) == 0) .OR. (ISS%hmask(i-1,j) == 2) ) then
! left face of the cell is at a stress boundary
! the depth-integrated longitudinal stress is equal to the difference of depth-integrated
! pressure on either side of the face
Expand All @@ -1981,19 +1988,19 @@ subroutine calc_shelf_driving_stress(CS, ISS, G, US, taudx, taudy, OD)
taudx(I-1,J) = taudx(I-1,J) - .5 * dyh * neumann_val
endif

if ((CS%u_face_mask(I,j) == 2) .OR. (ISS%hmask(i+1,j) == 0) .OR. (ISS%hmask(i+1,j) == 2) ) then
if ((CS%u_face_mask_bdry(I,j) == 2) .OR. (ISS%hmask(i+1,j) == 0) .OR. (ISS%hmask(i+1,j) == 2) ) then
! east face of the cell is at a stress boundary
taudx(I,J-1) = taudx(I,J-1) + .5 * dyh * neumann_val
taudx(I,J) = taudx(I,J) + .5 * dyh * neumann_val
endif

if ((CS%v_face_mask(i,J-1) == 2) .OR. (ISS%hmask(i,j-1) == 0) .OR. (ISS%hmask(i,j-1) == 2) ) then
if ((CS%v_face_mask_bdry(i,J-1) == 2) .OR. (ISS%hmask(i,j-1) == 0) .OR. (ISS%hmask(i,j-1) == 2) ) then
! south face of the cell is at a stress boundary
taudy(I-1,J-1) = taudy(I-1,J-1) - .5 * dxh * neumann_val
taudy(I,J-1) = taudy(I,J-1) - .5 * dxh * neumann_val
endif

if ((CS%v_face_mask(i,J) == 2) .OR. (ISS%hmask(i,j+1) == 0) .OR. (ISS%hmask(i,j+1) == 2) ) then
if ((CS%v_face_mask_bdry(i,J) == 2) .OR. (ISS%hmask(i,j+1) == 0) .OR. (ISS%hmask(i,j+1) == 2) ) then
! north face of the cell is at a stress boundary
taudy(I-1,J) = taudy(I-1,J) + .5 * dxh * neumann_val
taudy(I,J) = taudy(I,J) + .5 * dxh * neumann_val
Expand Down Expand Up @@ -2560,7 +2567,7 @@ subroutine apply_boundary_values(CS, ISS, G, US, time, Phisub, H_node, ice_visc,
end subroutine apply_boundary_values


!> Update depth integrated viscosity, based on horizontal strain rates, and also update the
!> Update depth integrated viscosity, based on horizontal strain rates
subroutine calc_shelf_visc(CS, ISS, G, US, u_shlf, v_shlf)
type(ice_shelf_dyn_CS), intent(inout) :: CS !< A pointer to the ice shelf control structure
type(ice_shelf_state), intent(in) :: ISS !< A structure with elements that describe
Expand All @@ -2575,7 +2582,6 @@ subroutine calc_shelf_visc(CS, ISS, G, US, u_shlf, v_shlf)
! quadrature points surrounding the cell vertices [L-1 ~> m-1].
! update DEPTH_INTEGRATED viscosity, based on horizontal strain rates - this is for bilinear FEM solve

! also this subroutine updates the nonlinear part of the basal traction

! this may be subject to change later... to make it "hybrid"
! real, dimension(SZDIB_(G),SZDJB_(G)) :: eII, ux, uy, vx, vy
Expand Down Expand Up @@ -2609,8 +2615,8 @@ subroutine calc_shelf_visc(CS, ISS, G, US, u_shlf, v_shlf)
do j=jsc,jec ; do i=isc,iec

if ((ISS%hmask(i,j) == 1) .OR. (ISS%hmask(i,j) == 3)) then
Visc_coef = US%kg_m2s_to_RZ_T*US%m_to_L*US%Z_to_L*(CS%AGlen_visc(i,j))**(-1./CS%n_glen)

Visc_coef = ( (US%RL2_T2_to_Pa)**(-CS%n_glen)*US%T_to_s )**(-1./CS%n_glen) * (CS%AGlen_visc(i,j))**(-1./CS%n_glen)
! Units of Aglen_visc [Pa-3 s-1]
do iq=1,2 ; do jq=1,2

ux = ( (u_shlf(I-1,J-1) * Phi(1,2*(jq-1)+iq,i,j) + &
Expand All @@ -2633,14 +2639,23 @@ subroutine calc_shelf_visc(CS, ISS, G, US, u_shlf, v_shlf)
(v_shlf(I-1,J) * Phi(6,2*(jq-1)+iq,i,j) + &
v_shlf(I,J-1) * Phi(4,2*(jq-1)+iq,i,j)) )
enddo ; enddo
! 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
if (trim(CS%ice_viscosity_compute)=="CONSTANT") then
CS%ice_visc(i,j) =1e15 * US%kg_m3_to_R*US%m_to_L*US%m_s_to_L_T * (G%areaT(i,j) * ISS%h_shelf(i,j))
! constant viscocity for debugging
elseif (trim(CS%ice_viscosity_compute)=="MODEL") then
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))
elseif(trim(CS%ice_viscosity_compute)=="OBS") then
if (CS%AGlen_visc(i,j) >0) CS%ice_visc(i,j) =CS%AGlen_visc(i,j)*(G%areaT(i,j) * ISS%h_shelf(i,j))
! Here CS%Aglen_visc(i,j) is the ice viscocity [Pa s-1] computed from obs and read from a file
endif
endif
enddo ; enddo
deallocate(Phi)
end subroutine calc_shelf_visc


!> Update basal shear
subroutine calc_shelf_taub(CS, ISS, G, US, u_shlf, v_shlf)
type(ice_shelf_dyn_CS), intent(inout) :: CS !< A pointer to the ice shelf control structure
type(ice_shelf_state), intent(in) :: ISS !< A structure with elements that describe
Expand Down
24 changes: 12 additions & 12 deletions src/ice_shelf/MOM_ice_shelf_initialize.F90
Original file line number Diff line number Diff line change
Expand Up @@ -395,14 +395,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,&
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(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)), &
intent(inout) :: v_shelf !< The meridional ice shelf velocity [L T-1 ~> m s-1].

real, dimension(SZDI_(G),SZDJ_(G)), &
intent(inout) :: float_cond !< An array indicating where the ice
!! shelf is floating: 0 if floating, 1 if not. [nondim]
Expand Down Expand Up @@ -450,14 +451,12 @@ subroutine initialize_ice_flow_from_file(bed_elev,u_shelf, v_shelf,float_cond,&

floatfr_varname = "float_frac"

!### I think that the following two lines should have ..., scale=US%m_s_to_L_T
call MOM_read_data(filename, trim(ushelf_varname), u_shelf, G%Domain, position=CORNER, scale=1.0)
call MOM_read_data(filename, trim(vshelf_varname), v_shelf, G%Domain, position=CORNER, scale=1.0)
! call MOM_read_data(filename, trim(ice_visc_varname), ice_visc, G%Domain,position=CORNER,scale=1.0)
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.)
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 @@ -543,11 +542,12 @@ subroutine initialize_ice_shelf_boundary_from_file(u_face_mask_bdry, v_face_mask
" initialize_ice_shelf_velocity_from_file: Unable to open "//trim(filename))


call MOM_read_data(filename, trim(ufcmskbdry_varname), u_face_mask_bdry, G%Domain, position=CORNER, scale=1.0)
call MOM_read_data(filename, trim(vfcmskbdry_varname), v_face_mask_bdry, G%Domain, position=CORNER, scale=1.0)
!### I think that the following two lines should have ..., scale=US%m_s_to_L_T
call MOM_read_data(filename, trim(ubdryv_varname), u_bdry_val, G%Domain, position=CORNER, scale=1.0)
call MOM_read_data(filename, trim(vbdryv_varname), v_bdry_val, G%Domain, position=CORNER, scale=1.)
call MOM_read_data(filename, trim(ufcmskbdry_varname), u_face_mask_bdry, G%Domain, position=CORNER, &
scale=US%m_s_to_L_T)
call MOM_read_data(filename, trim(vfcmskbdry_varname), v_face_mask_bdry, G%Domain, position=CORNER, &
scale=US%m_s_to_L_T)
call MOM_read_data(filename, trim(ubdryv_varname), u_bdry_val, G%Domain, position=CORNER, scale=US%m_s_to_L_T)
call MOM_read_data(filename, trim(vbdryv_varname), v_bdry_val, G%Domain, position=CORNER, scale=US%m_s_to_L_T)
call MOM_read_data(filename, trim(umask_varname), umask, G%Domain, position=CORNER, scale=1.)
call MOM_read_data(filename, trim(vmask_varname), vmask, G%Domain, position=CORNER, scale=1.)
filename = trim(inputdir)//trim(icethick_file)
Expand Down Expand Up @@ -615,7 +615,7 @@ subroutine initialize_ice_C_basal_friction(C_basal_friction, G, US, PF)
end subroutine


!> Initialize ice basal friction
!> Initialize ice-stiffness parameter
subroutine initialize_ice_AGlen(AGlen, G, US, PF)
type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure
real, dimension(SZDI_(G),SZDJ_(G)), &
Expand Down

0 comments on commit f92a4ab

Please sign in to comment.