Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Ice dynamics #1479

Merged
merged 13 commits into from
Sep 1, 2021
Merged

Ice dynamics #1479

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
13 commits
Select commit Hold shift + click to select a range
7eb92f6
Merge pull request #2 from NOAA-GFDL/dev/gfdl
OlgaSergienko Apr 15, 2021
13f5e12
Merge https://github.com/OlgaSergienko/MOM6 into ice_dynamics
OlgaSergienko Jun 17, 2021
457c499
increased a number of iterations in solver for non-linear viscosity
OlgaSergienko Jun 17, 2021
0267d08
Merge branch 'ice_dynamics' of https://github.com/OlgaSergienko/MOM6 …
OlgaSergienko Jun 17, 2021
347e147
Merge branch 'NOAA-GFDL:dev/gfdl' into dev/gfdl
OlgaSergienko Jun 17, 2021
016beb6
Merge branch 'NOAA-GFDL:dev/gfdl' into dev/gfdl
OlgaSergienko Jul 26, 2021
8ef9596
Merge https://github.com/OlgaSergienko/MOM6 into ice_dynamics
OlgaSergienko Jul 26, 2021
3ce2efe
changes to ice_velocity_mask_update, initialization of ice velocity a…
OlgaSergienko Aug 25, 2021
a5caf2e
removed tab from MOM_ice_shelf_initialize.F90
OlgaSergienko Aug 25, 2021
442ae8f
changed the block idetentation for consistency with MOM6 code style
OlgaSergienko Aug 31, 2021
e2c565c
changed identation MOM_ice_shelf_initialize.F90 to follow MOM6 conven…
OlgaSergienko Aug 31, 2021
f46c4bc
Changed identation in MOM6_ice_shelf_initialize.F90 to follow MOM6 co…
OlgaSergienko Sep 1, 2021
d234258
Merge branch 'dev/gfdl' into ice_dynamics
marshallward Sep 1, 2021
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
159 changes: 79 additions & 80 deletions src/ice_shelf/MOM_ice_shelf_dynamics.F90
View file
Open in desktop
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
Hallberg-NOAA marked this conversation as resolved.
Show resolved Hide resolved
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
Open in desktop
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
Hallberg-NOAA marked this conversation as resolved.
Show resolved Hide resolved
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