From 1fcc4e77717818e0826a3347724cc1e259315d8d Mon Sep 17 00:00:00 2001 From: alex-huth Date: Wed, 5 Jun 2024 14:09:27 -0400 Subject: [PATCH 1/6] Updated with dev/gfdl. Then, added parameters for min ice thickness, min basal traction, max surf slope, and min ice viscosity to use for ice dynamics --- src/ice_shelf/MOM_ice_shelf_dynamics.F90 | 80 ++++++++++++++++-------- 1 file changed, 55 insertions(+), 25 deletions(-) diff --git a/src/ice_shelf/MOM_ice_shelf_dynamics.F90 b/src/ice_shelf/MOM_ice_shelf_dynamics.F90 index bee5cf11aa..a1ed8512fa 100644 --- a/src/ice_shelf/MOM_ice_shelf_dynamics.F90 +++ b/src/ice_shelf/MOM_ice_shelf_dynamics.F90 @@ -107,7 +107,7 @@ module MOM_ice_shelf_dynamics !! of "linearized" basal stress (Pa) [R L3 T-1 ~> kg s-1] !! 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 s-1)^(n_basal_fric) + !! units= Pa (s m-1)^(n_basal_fric) real, pointer, dimension(:,:) :: OD_rt => NULL() !< A running total for calculating OD_av [Z ~> m]. real, pointer, dimension(:,:) :: ground_frac_rt => NULL() !< A running total for calculating ground_frac. real, pointer, dimension(:,:) :: OD_av => NULL() !< The time average open ocean depth [Z ~> m]. @@ -164,6 +164,11 @@ module MOM_ice_shelf_dynamics real :: CFL_factor !< A factor used to limit subcycled advective timestep in uncoupled runs !! i.e. dt <= CFL_factor * min(dx / u) [nondim] + real :: min_h_shelf !< The minimum ice thickness used during ice dynamics [L ~> m]. + real :: min_basal_traction !< The minimum basal traction for grounded ice (Pa m-1 s) [R L T-1 ~> kg m-2 s-1] + real :: max_surface_slope !< The maximum allowed ice-sheet surface slope (to ignore, set to zero) [nondim] + real :: min_ice_visc !< The minimum allowed Glen's law ice viscosity (Pa s), in [R L2 T-1 ~> kg m-1 s-1]. + real :: n_glen !< Nonlinearity exponent in Glen's Law [nondim] real :: eps_glen_min !< Min. strain rate to avoid infinite Glen's law viscosity, [T-1 ~> s-1]. real :: n_basal_fric !< Exponent in sliding law tau_b = C u^(m_slide) [nondim] @@ -343,7 +348,7 @@ subroutine register_ice_shelf_dyn_restarts(G, US, param_file, CS, restart_CS) allocate(CS%ice_visc(isd:ied,jsd:jed,CS%visc_qps), 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) ! [R L3 T-1 ~> kg s-1] - allocate(CS%C_basal_friction(isd:ied,jsd:jed), source=5.0e10) ! [Pa (m-1 s)^n_sliding] + allocate(CS%C_basal_friction(isd:ied,jsd:jed), source=5.0e10) ! [Pa (s m-1)^n_sliding] allocate(CS%OD_av(isd:ied,jsd:jed), source=0.0) allocate(CS%ground_frac(isd:ied,jsd:jed), source=0.0) allocate(CS%taudx_shelf(IsdB:IedB,JsdB:JedB), source=0.0) @@ -378,7 +383,7 @@ subroutine register_ice_shelf_dyn_restarts(G, US, param_file, CS, restart_CS) 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, "C_basal_friction", .true., restart_CS, & - "basal sliding coefficients", "Pa (m s-1)^n_sliding") + "basal sliding coefficients", "Pa (s m-1)^n_sliding") 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_val", .false., restart_CS, & @@ -491,6 +496,19 @@ subroutine initialize_ice_shelf_dyn(param_file, Time, ISS, CS, G, US, diag, new_ "The gravitational acceleration of the Earth.", & units="m s-2", default=9.80, scale=US%m_s_to_L_T**2*US%Z_to_m) + call get_param(param_file, mdl, "MIN_H_SHELF", CS%min_h_shelf, & + "min. ice thickness used during ice dynamics", & + units="m", default=0.,scale=US%m_to_L) + call get_param(param_file, mdl, "MIN_BASAL_TRACTION", CS%min_basal_traction, & + "min. allowed basal traction", & + units="Pa m-1 yr", default=0., scale=365.0*86400.0*US%Pa_to_RLZ_T2*US%L_T_to_m_s) + call get_param(param_file, mdl, "MAX_SURFACE_SLOPE", CS%max_surface_slope, & + "max. allowed ice-sheet surface slope. To ignore, set to zero.", & + units="none", default=0., scale=US%m_to_Z/US%m_to_L) + call get_param(param_file, mdl, "MIN_ICE_VISC", CS%min_ice_visc, & + "min. allowed Glen's law ice viscosity", & + units="Pa s", default=0., scale=US%Pa_to_RL2_T2*US%s_to_T) + call get_param(param_file, mdl, "GLEN_EXPONENT", CS%n_glen, & "nonlinearity exponent in Glen's Law", & units="none", default=3.) @@ -837,7 +855,7 @@ subroutine initialize_diagnostic_fields(CS, ISS, G, US, Time) do j=jsd,jed do i=isd,ied - OD = CS%bed_elev(i,j) - rhoi_rhow * ISS%h_shelf(i,j) + OD = CS%bed_elev(i,j) - rhoi_rhow * max(ISS%h_shelf(i,j),CS%min_h_shelf) if (OD >= 0) then ! ice thickness does not take up whole ocean column -> floating CS%OD_av(i,j) = OD @@ -1328,7 +1346,7 @@ subroutine ice_shelf_solve_outer(CS, ISS, G, US, u_shlf, v_shlf, taudx, taudy, i if (.not. CS%GL_couple) then do j=G%jsc,G%jec ; do i=G%isc,G%iec - if (rhoi_rhow * ISS%h_shelf(i,j) - CS%bed_elev(i,j) > 0) then + if (rhoi_rhow * max(ISS%h_shelf(i,j),CS%min_h_shelf) - CS%bed_elev(i,j) > 0) then CS%ground_frac(i,j) = 1.0 CS%OD_av(i,j) =0.0 endif @@ -1346,7 +1364,7 @@ subroutine ice_shelf_solve_outer(CS, ISS, G, US, u_shlf, v_shlf, taudx, taudy, i if (CS%GL_regularize) then - call interpolate_H_to_B(G, ISS%h_shelf, ISS%hmask, H_node) + call interpolate_H_to_B(G, ISS%h_shelf, ISS%hmask, H_node, CS%min_h_shelf) do j=G%jsc,G%jec ; do i=G%isc,G%iec nodefloat = 0 @@ -2263,7 +2281,7 @@ subroutine calc_shelf_driving_stress(CS, ISS, G, US, taudx, taudy, OD) real :: neumann_val ! [R Z L2 T-2 ~> kg s-2] real :: dxh, dyh,Dx,Dy ! Local grid spacing [L ~> m] real :: grav ! The gravitational acceleration [L2 Z-1 T-2 ~> m s-2] - + real :: scale ! Scaling factor used to ensure surface slope magnitude does not exceed CS%max_surface_slope integer :: i, j, iscq, iecq, jscq, jecq, isd, jsd, ied, jed, is, js, iegq, jegq integer :: giec, gjec, gisc, gjsc, cnt, isc, jsc, iec, jec integer :: i_off, j_off @@ -2289,17 +2307,17 @@ subroutine calc_shelf_driving_stress(CS, ISS, G, US, taudx, taudy, OD) if (CS%GL_couple) then do j=jsc-G%domain%njhalo,jec+G%domain%njhalo do i=isc-G%domain%nihalo,iec+G%domain%nihalo - S(i,j) = -CS%bed_elev(i,j) + (OD(i,j) + ISS%h_shelf(i,j)) + S(i,j) = -CS%bed_elev(i,j) + (OD(i,j) + max(ISS%h_shelf(i,j),CS%min_h_shelf)) enddo enddo else ! check whether the ice is floating or grounded do j=jsc-G%domain%njhalo,jec+G%domain%njhalo do i=isc-G%domain%nihalo,iec+G%domain%nihalo - if (rhoi_rhow * ISS%h_shelf(i,j) - CS%bed_elev(i,j) <= 0) then - S(i,j) = (1 - rhoi_rhow)*ISS%h_shelf(i,j) + if (rhoi_rhow * max(ISS%h_shelf(i,j),CS%min_h_shelf) - CS%bed_elev(i,j) <= 0) then + S(i,j) = (1 - rhoi_rhow)*max(ISS%h_shelf(i,j),CS%min_h_shelf) else - S(i,j) = ISS%h_shelf(i,j)-CS%bed_elev(i,j) + S(i,j) = max(ISS%h_shelf(i,j),CS%min_h_shelf)-CS%bed_elev(i,j) endif enddo enddo @@ -2393,14 +2411,19 @@ subroutine calc_shelf_driving_stress(CS, ISS, G, US, taudx, taudy, OD) endif endif - sx_e(i,j) = (-.25 * G%areaT(i,j)) * ((rho * grav) * (ISS%h_shelf(i,j) * sx)) - sy_e(i,j) = (-.25 * G%areaT(i,j)) * ((rho * grav) * (ISS%h_shelf(i,j) * sy)) + if (CS%max_surface_slope>0) then + scale = min(CS%max_surface_slope/sqrt(sx**2+sy**2),1.0) + sx = scale*sx; sy = scale*sy + endif + + sx_e(i,j) = (-.25 * G%areaT(i,j)) * ((rho * grav) * (max(ISS%h_shelf(i,j),CS%min_h_shelf) * sx)) + sy_e(i,j) = (-.25 * G%areaT(i,j)) * ((rho * grav) * (max(ISS%h_shelf(i,j),CS%min_h_shelf) * sy)) !Stress (Neumann) boundary conditions 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) * (rho * max(ISS%h_shelf(i,j),CS%min_h_shelf)**2 - rhow * CS%bed_elev(i,j)**2)) else - neumann_val = (.5 * grav) * ((1-rho/rhow) * (rho * ISS%h_shelf(i,j)**2)) + neumann_val = (.5 * grav) * ((1-rho/rhow) * (rho * max(ISS%h_shelf(i,j),CS%min_h_shelf)**2)) endif if ((CS%u_face_mask_bdry(I-1,j) == 2) .OR. & ((ISS%hmask(i-1,j) == 0 .OR. ISS%hmask(i-1,j) == 2) .AND. (CS%reentrant_x .OR. (i+i_off /= gisc)))) then @@ -2996,10 +3019,14 @@ subroutine calc_shelf_visc(CS, ISS, G, US, u_shlf, v_shlf) if ((ISS%hmask(i,j) == 1) .OR. (ISS%hmask(i,j) == 3)) then if (trim(CS%ice_viscosity_compute) == "CONSTANT") then - CS%ice_visc(i,j,1) = 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)) + CS%ice_visc(i,j,1) = 1e15 * (US%kg_m3_to_R*US%m_to_L*US%m_s_to_L_T) * & + (G%areaT(i,j) * max(ISS%h_shelf(i,j),CS%min_h_shelf)) ! constant viscocity for debugging elseif (trim(CS%ice_viscosity_compute) == "OBS") then - if (CS%AGlen_visc(i,j) >0) CS%ice_visc(i,j,1) = CS%AGlen_visc(i,j) * (G%areaT(i,j) * ISS%h_shelf(i,j)) + if (CS%AGlen_visc(i,j) >0) then + CS%ice_visc(i,j,1) = max(CS%AGlen_visc(i,j) * (G%areaT(i,j) * max(ISS%h_shelf(i,j),CS%min_h_shelf)),& + CS%min_ice_visc) + endif ! Here CS%Aglen_visc(i,j) is the ice viscosity [Pa s ~> R L2 T-1] computed from obs and read from a file elseif (model_qp1) then !calculate viscosity at 1 cell-centered quadrature point per cell @@ -3027,9 +3054,9 @@ subroutine calc_shelf_visc(CS, ISS, G, US, u_shlf, v_shlf) (v_shlf(I-1,J) * CS%PhiC(6,i,j) + & v_shlf(I,J-1) * CS%PhiC(4,i,j)) - CS%ice_visc(i,j,1) = 0.5 * Visc_coef * (G%areaT(i,j) * ISS%h_shelf(i,j)) * & + CS%ice_visc(i,j,1) = max(0.5 * Visc_coef * (G%areaT(i,j) * max(ISS%h_shelf(i,j),CS%min_h_shelf)) * & (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)) * & - (US%Pa_to_RL2_T2*US%s_to_T) + (US%Pa_to_RL2_T2*US%s_to_T),CS%min_ice_visc) elseif (model_qp4) then !calculate viscosity at 4 quadrature points per cell @@ -3057,9 +3084,9 @@ subroutine calc_shelf_visc(CS, ISS, G, US, u_shlf, v_shlf) (v_shlf(I,J-1) * CS%Phi(4,2*(jq-1)+iq,i,j) + & v_shlf(I-1,J) * CS%Phi(6,2*(jq-1)+iq,i,j)) - CS%ice_visc(i,j,2*(jq-1)+iq) = 0.5 * Visc_coef * (G%areaT(i,j) * ISS%h_shelf(i,j)) * & + CS%ice_visc(i,j,2*(jq-1)+iq) = max(0.5 * Visc_coef * (G%areaT(i,j) * max(ISS%h_shelf(i,j),CS%min_h_shelf)) * & (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)) * & - (US%Pa_to_RL2_T2*US%s_to_T) + (US%Pa_to_RL2_T2*US%s_to_T),CS%min_ice_visc) enddo; enddo endif endif @@ -3123,7 +3150,7 @@ subroutine calc_shelf_taub(CS, ISS, G, US, u_shlf, v_shlf) if (CS%CoulombFriction) then !Effective pressure Hf = max((CS%density_ocean_avg/CS%density_ice) * CS%bed_elev(i,j), 0.0) - fN = max(fN_scale*((CS%density_ice * CS%g_Earth) * (ISS%h_shelf(i,j) - Hf)),CS%CF_MinN) + fN = max(fN_scale*((CS%density_ice * CS%g_Earth) * (max(ISS%h_shelf(i,j),CS%min_h_shelf) - Hf)),CS%CF_MinN) fB = alpha * (CS%C_basal_friction(i,j) / (CS%CF_Max * fN))**(CS%CF_PostPeak/CS%n_basal_fric) CS%basal_traction(i,j) = ((G%areaT(i,j) * CS%C_basal_friction(i,j)) * & @@ -3134,6 +3161,8 @@ subroutine calc_shelf_taub(CS, ISS, G, US, u_shlf, v_shlf) CS%basal_traction(i,j) = ((G%areaT(i,j) * CS%C_basal_friction(i,j)) * (unorm**(CS%n_basal_fric-1))) * & (US%Pa_to_RLZ_T2*US%L_T_to_m_s) endif + + CS%basal_traction(i,j)=max(CS%basal_traction(i,j), CS%min_basal_traction * G%areaT(i,j)) endif enddo enddo @@ -3194,7 +3223,7 @@ subroutine update_OD_ffrac_uncoupled(CS, G, h_shelf) do j=jsd,jed do i=isd,ied - OD = CS%bed_elev(i,j) - rhoi_rhow * h_shelf(i,j) + OD = CS%bed_elev(i,j) - rhoi_rhow * max(h_shelf(i,j),CS%min_h_shelf) if (OD >= 0) then ! ice thickness does not take up whole ocean column -> floating CS%OD_av(i,j) = OD @@ -3640,7 +3669,7 @@ end subroutine update_velocity_masks !> Interpolate the ice shelf thickness from tracer point to nodal points, !! subject to a mask. -subroutine interpolate_H_to_B(G, h_shelf, hmask, H_node) +subroutine interpolate_H_to_B(G, h_shelf, hmask, H_node, min_h_shelf) type(ocean_grid_type), intent(inout) :: G !< The grid structure used by the ice shelf. real, dimension(SZDI_(G),SZDJ_(G)), & intent(in) :: h_shelf !< The ice shelf thickness at tracer points [Z ~> m]. @@ -3650,6 +3679,7 @@ subroutine interpolate_H_to_B(G, h_shelf, hmask, H_node) real, dimension(SZDIB_(G),SZDJB_(G)), & intent(inout) :: H_node !< The ice shelf thickness at nodal (corner) !! points [Z ~> m]. + real, intent(in) :: min_h_shelf !< The minimum ice thickness used during ice dynamics [L ~> m]. integer :: i, j, isc, iec, jsc, jec, num_h, k, l, ic, jc real :: h_arr(2,2) @@ -3666,7 +3696,7 @@ subroutine interpolate_H_to_B(G, h_shelf, hmask, H_node) num_h = 0 do l=1,2; jc=j-1+l; do k=1,2; ic=i-1+k if (hmask(ic,jc) == 1.0 .or. hmask(ic,jc) == 3.0) then - h_arr(k,l)=h_shelf(ic,jc) + h_arr(k,l)=max(h_shelf(ic,jc),min_h_shelf) num_h = num_h + 1 else h_arr(k,l)=0.0 From c9810fbabe0b333ff4baa025269c1fb524603a9f Mon Sep 17 00:00:00 2001 From: alex-huth Date: Thu, 6 Jun 2024 17:00:23 -0400 Subject: [PATCH 2/6] noted how the MIN_BASAL_TRACTION parameter input is in units [Pa m-1 yr], but converts automatically to [Pa m-1 s] in the code --- src/ice_shelf/MOM_ice_shelf_dynamics.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/ice_shelf/MOM_ice_shelf_dynamics.F90 b/src/ice_shelf/MOM_ice_shelf_dynamics.F90 index a1ed8512fa..c9a3baa8e2 100644 --- a/src/ice_shelf/MOM_ice_shelf_dynamics.F90 +++ b/src/ice_shelf/MOM_ice_shelf_dynamics.F90 @@ -500,7 +500,7 @@ subroutine initialize_ice_shelf_dyn(param_file, Time, ISS, CS, G, US, diag, new_ "min. ice thickness used during ice dynamics", & units="m", default=0.,scale=US%m_to_L) call get_param(param_file, mdl, "MIN_BASAL_TRACTION", CS%min_basal_traction, & - "min. allowed basal traction", & + "min. allowed basal traction. Input is in [Pa m-1 yr], but is converted when read in to [Pa m-1 s]", & units="Pa m-1 yr", default=0., scale=365.0*86400.0*US%Pa_to_RLZ_T2*US%L_T_to_m_s) call get_param(param_file, mdl, "MAX_SURFACE_SLOPE", CS%max_surface_slope, & "max. allowed ice-sheet surface slope. To ignore, set to zero.", & From 066c5e236e79837816b129c7b3fc488745f77608 Mon Sep 17 00:00:00 2001 From: alex-huth Date: Wed, 26 Jun 2024 17:52:38 -0400 Subject: [PATCH 3/6] Added separate scalar ice-shelf diagnostics for Antarctica and Greenland --- src/ice_shelf/MOM_ice_shelf.F90 | 498 +++++++++++++++++------ src/ice_shelf/MOM_ice_shelf_dynamics.F90 | 28 +- 2 files changed, 398 insertions(+), 128 deletions(-) diff --git a/src/ice_shelf/MOM_ice_shelf.F90 b/src/ice_shelf/MOM_ice_shelf.F90 index 4f811cac87..18e3669f98 100644 --- a/src/ice_shelf/MOM_ice_shelf.F90 +++ b/src/ice_shelf/MOM_ice_shelf.F90 @@ -181,6 +181,8 @@ module MOM_ice_shelf !! fluxes. It will avoid large increase in sea level. logical :: constant_sea_level_misomip !< If true, constant_sea_level fluxes are applied only over !! the surface sponge cells from the ISOMIP/MISOMIP configuration + logical :: smb_diag !< If true, calculate diagnostics related to surface mass balance + logical :: bmb_diag !< If true, calculate diagnostics related to basal mass balance real :: min_ocean_mass_float !< The minimum ocean mass per unit area before the ice !! shelf is considered to float when constant_sea_level !! is used [R Z ~> kg m-2] @@ -213,7 +215,17 @@ module MOM_ice_shelf id_bdott_melt = -1, id_bdott_accum = -1, id_bdott = -1, & id_dvafdt = -1, id_g_adot = -1, id_f_adot = -1, id_adot = -1, & id_bdot_melt = -1, id_bdot_accum = -1, id_bdot = -1, & - id_t_area = -1, id_g_area = -1, id_f_area = -1 + id_t_area = -1, id_g_area = -1, id_f_area = -1, & + id_Ant_vaf = -1, id_Ant_g_adott = -1, id_Ant_f_adott = -1, id_Ant_adott = -1, & + id_Ant_bdott_melt = -1, id_Ant_bdott_accum = -1, id_Ant_bdott = -1, & + id_Ant_dvafdt = -1, id_Ant_g_adot = -1, id_Ant_f_adot = -1, id_Ant_adot = -1, & + id_Ant_bdot_melt = -1, id_Ant_bdot_accum = -1, id_Ant_bdot = -1, & + id_Ant_t_area = -1, id_Ant_g_area = -1, id_Ant_f_area = -1, & + id_Gr_vaf = -1, id_Gr_g_adott = -1, id_Gr_f_adott = -1, id_Gr_adott = -1, & + id_Gr_bdott_melt = -1, id_Gr_bdott_accum = -1, id_Gr_bdott = -1, & + id_Gr_dvafdt = -1, id_Gr_g_adot = -1, id_Gr_f_adot = -1, id_Gr_adot = -1, & + id_Gr_bdot_melt = -1, id_Gr_bdot_accum = -1, id_Gr_bdot = -1, & + id_Gr_t_area = -1, id_Gr_g_area = -1, id_Gr_f_area = -1 !>@} type(external_field) :: mass_handle @@ -270,12 +282,10 @@ subroutine shelf_calc_flux(sfc_state_in, fluxes_in, Time, time_step_in, CS) p_int !< The pressure at the ice-ocean interface [R L2 T-2 ~> Pa]. real, dimension(SZI_(CS%grid),SZJ_(CS%grid)) :: & - exch_vel_t, & !< Sub-shelf thermal exchange velocity [Z T-1 ~> m s-1] - exch_vel_s, & !< Sub-shelf salt exchange velocity [Z T-1 ~> m s-1] - tmp, & !< Temporary field used when calculating diagnostics [various] - dh_bdott, & !< Basal melt/accumulation over a time step, used for diagnostics [Z ~> m] - dh_adott !< Surface melt/accumulation over a time step, used for diagnostics [Z ~> m] - + exch_vel_t, & !< Sub-shelf thermal exchange velocity [Z T-1 ~> m s-1] + exch_vel_s, & !< Sub-shelf salt exchange velocity [Z T-1 ~> m s-1] + dh_bdott, & !< Basal melt/accumulation over a time step, used for diagnostics [Z ~> m] + dh_adott !< Surface melt/accumulation over a time step, used for diagnostics [Z ~> m] real, dimension(SZDI_(CS%grid),SZDJ_(CS%grid)) :: & mass_flux !< Total mass flux of freshwater across the ice-ocean interface. [R Z L2 T-1 ~> kg s-1] real, dimension(SZDI_(CS%grid),SZDJ_(CS%grid)) :: & @@ -343,9 +353,8 @@ subroutine shelf_calc_flux(sfc_state_in, fluxes_in, Time, time_step_in, CS) character(len=160) :: mesg ! The text of an error message integer, dimension(2) :: EOSdom ! The i-computational domain for the equation of state integer :: i, j, is, ie, js, je, ied, jed, it1, it3 - real :: vaf0, vaf ! The previous and current volume above floatation [m3] - logical :: smb_diag=.false., bmb_diag=.false. ! Flags to calculate diagnostics related to surface/basal mass balance - real :: val ! Temporary value when calculating scalar diagnostics [various] + real :: vaf0, vaf0_A, vaf0_G !The previous volumes above floatation [m3] + !for all ice sheets, Antarctica only, or Greenland only [m3] if (.not. associated(CS)) call MOM_error(FATAL, "shelf_calc_flux: "// & "initialize_ice_shelf must be called before shelf_calc_flux.") @@ -356,13 +365,14 @@ subroutine shelf_calc_flux(sfc_state_in, fluxes_in, Time, time_step_in, CS) time_step = time_step_in Itime_step = 1./time_step - if (CS%id_adott>0 .or. CS%id_g_adott>0 .or. CS%id_f_adott>0 .or. & - CS%id_adot >0 .or. CS%id_g_adot >0 .or. CS%id_f_adot >0 ) smb_diag=.true. - if (CS%id_bdott>0 .or. CS%id_bdott_melt>0 .or. CS%id_bdott_accum>0 .or. & - CS%id_bdot >0 .or. CS%id_bdot_melt >0 .or. CS%id_bdot_accum >0) bmb_diag=.true. + dh_adott(:,:)=0.0; dh_bdott(:,:)=0.0 - if (CS%active_shelf_dynamics .and. CS%id_dvafdt > 0) & !calculate previous volume above floatation - call volume_above_floatation(CS%dCS, G, ISS, vaf0) + if (CS%active_shelf_dynamics) then + !calculate previous volumes above floatation + if (CS%id_dvafdt > 0) call volume_above_floatation(CS%dCS, G, ISS, vaf0) !all ice sheet + if (CS%id_Ant_dvafdt > 0) call volume_above_floatation(CS%dCS, G, ISS, vaf0_A, hemisphere=0) !Antarctica only + if (CS%id_Gr_dvafdt > 0) call volume_above_floatation(CS%dCS, G, ISS, vaf0_G, hemisphere=1) !Greenland only + endif is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; ied = G%ied ; jed = G%jed if (CS%data_override_shelf_fluxes .and. CS%active_shelf_dynamics) then @@ -766,9 +776,9 @@ subroutine shelf_calc_flux(sfc_state_in, fluxes_in, Time, time_step_in, CS) ! Melting has been computed, now is time to update thickness and mass if ( CS%override_shelf_movement .and. (.not.CS%mass_from_file)) then - if (bmb_diag) dh_bdott(is:ie,js:je) = ISS%h_shelf(is:ie,js:je) + if (CS%bmb_diag) dh_bdott(is:ie,js:je) = ISS%h_shelf(is:ie,js:je) call change_thickness_using_melt(ISS, G, US, time_step, fluxes, CS%density_ice, CS%debug) - if (bmb_diag) dh_bdott(is:ie,js:je) = ISS%h_shelf(is:ie,js:je) - dh_bdott(is:ie,js:je) + if (CS%bmb_diag) dh_bdott(is:ie,js:je) = ISS%h_shelf(is:ie,js:je) - dh_bdott(is:ie,js:je) if (CS%debug) then call hchksum(ISS%h_shelf, "h_shelf after change thickness using melt", G%HI, haloshift=0, unscale=US%Z_to_m) @@ -782,9 +792,9 @@ subroutine shelf_calc_flux(sfc_state_in, fluxes_in, Time, time_step_in, CS) ISS%dhdt_shelf(:,:) = ISS%h_shelf(:,:) - if (bmb_diag) dh_bdott(is:ie,js:je) = ISS%h_shelf(is:ie,js:je) + if (CS%bmb_diag) dh_bdott(is:ie,js:je) = ISS%h_shelf(is:ie,js:je) call change_thickness_using_melt(ISS, G, US, time_step, fluxes, CS%density_ice, CS%debug) - if (bmb_diag) dh_bdott(is:ie,js:je) = ISS%h_shelf(is:ie,js:je) - dh_bdott(is:ie,js:je) + if (CS%bmb_diag) dh_bdott(is:ie,js:je) = ISS%h_shelf(is:ie,js:je) - dh_bdott(is:ie,js:je) if (CS%debug) then call hchksum(ISS%h_shelf, "h_shelf after change thickness using melt", G%HI, haloshift=0, unscale=US%Z_to_m) @@ -792,9 +802,9 @@ subroutine shelf_calc_flux(sfc_state_in, fluxes_in, Time, time_step_in, CS) unscale=US%RZ_to_kg_m2) endif - if (smb_diag) dh_adott(is:ie,js:je) = ISS%h_shelf(is:ie,js:je) + if (CS%smb_diag) dh_adott(is:ie,js:je) = ISS%h_shelf(is:ie,js:je) call change_thickness_using_precip(CS, ISS, G, US, fluxes, time_step, Time) - if (smb_diag) dh_adott(is:ie,js:je) = ISS%h_shelf(is:ie,js:je) - dh_adott(is:ie,js:je) + if (CS%smb_diag) dh_adott(is:ie,js:je) = ISS%h_shelf(is:ie,js:je) - dh_adott(is:ie,js:je) if (CS%debug) then call hchksum(ISS%h_shelf, "h_shelf after change thickness using surf acc", G%HI, haloshift=0, unscale=US%Z_to_m) @@ -846,69 +856,7 @@ subroutine shelf_calc_flux(sfc_state_in, fluxes_in, Time, time_step_in, CS) if (CS%id_h_shelf > 0) call post_data(CS%id_h_shelf, ISS%h_shelf, CS%diag) if (CS%id_dhdt_shelf > 0) call post_data(CS%id_dhdt_shelf, ISS%dhdt_shelf, CS%diag) if (CS%id_h_mask > 0) call post_data(CS%id_h_mask,ISS%hmask,CS%diag) - !scalars - if (CS%active_shelf_dynamics) then - if (CS%id_vaf > 0 .or. CS%id_dvafdt > 0) & !calculate current volume above floatation (vaf) - call volume_above_floatation(CS%dCS, G, ISS, vaf) - if (CS%id_vaf > 0) call post_scalar_data(CS%id_vaf ,vaf ,CS%diag) !current vaf - if (CS%id_dvafdt > 0) call post_scalar_data(CS%id_dvafdt,(vaf-vaf0)*Itime_step,CS%diag) !d(vaf)/dt - if (CS%id_adott > 0 .or. CS%id_adot > 0) then !surface accumulation - surface melt - call integrate_over_ice_sheet_area(G, ISS, dh_adott, US%Z_to_m, val) - if (CS%id_adott > 0) call post_scalar_data(CS%id_adott,val ,CS%diag) - if (CS%id_adot > 0) call post_scalar_data(CS%id_adot ,val*Itime_step,CS%diag) - endif - if (CS%id_g_adott > 0 .or. CS%id_g_adot > 0) then !grounded only: surface accumulation - surface melt - call masked_var_grounded(G,CS%dCS,dh_adott,tmp) - call integrate_over_ice_sheet_area(G, ISS, tmp, US%Z_to_m, val) - if (CS%id_g_adott > 0) call post_scalar_data(CS%id_g_adott,val ,CS%diag) - if (CS%id_g_adot > 0) call post_scalar_data(CS%id_g_adot ,val*Itime_step,CS%diag) - endif - if (CS%id_f_adott > 0 .or. CS%id_f_adot > 0) then !floating only: surface accumulation - surface melt - call masked_var_grounded(G,CS%dCS,dh_adott,tmp) - tmp(:,:) = dh_adott(:,:) - tmp(:,:) - call integrate_over_ice_sheet_area(G, ISS, tmp, US%Z_to_m, val) - if (CS%id_f_adott > 0) call post_scalar_data(CS%id_f_adott,val ,CS%diag) - if (CS%id_f_adot > 0) call post_scalar_data(CS%id_f_adot ,val*Itime_step,CS%diag) - endif - endif - if (CS%id_bdott > 0 .or. CS%id_bdot > 0) then !bottom accumulation - bottom melt - call integrate_over_ice_sheet_area(G, ISS, dh_bdott, US%Z_to_m, val) - if (CS%id_bdott > 0) call post_scalar_data(CS%id_bdott,val ,CS%diag) - if (CS%id_bdot > 0) call post_scalar_data(CS%id_bdot ,val*Itime_step,CS%diag) - endif - if (CS%id_bdott_melt > 0 .or. CS%id_bdot_melt > 0) then !bottom melt - tmp(:,:)=0.0 - do j=js,je ; do i=is,ie - if (dh_bdott(i,j) < 0) tmp(i,j) = -dh_bdott(i,j) - enddo; enddo - call integrate_over_ice_sheet_area(G, ISS, tmp, US%Z_to_m, val) - if (CS%id_bdott_melt > 0) call post_scalar_data(CS%id_bdott_melt,val ,CS%diag) - if (CS%id_bdot_melt > 0) call post_scalar_data(CS%id_bdot_melt ,val*Itime_step,CS%diag) - endif - if (CS%id_bdott_accum > 0 .or. CS%id_bdot_accum > 0) then !bottom accumulation - tmp(:,:)=0.0 - do j=js,je ; do i=is,ie - if (dh_bdott(i,j) > 0) tmp(i,j) = dh_bdott(i,j) - enddo; enddo - call integrate_over_ice_sheet_area(G, ISS, tmp, US%Z_to_m, val) - if (CS%id_bdott_accum > 0) call post_scalar_data(CS%id_bdott_accum,val ,CS%diag) - if (CS%id_bdot_accum > 0) call post_scalar_data(CS%id_bdot_accum ,val*Itime_step,CS%diag) - endif - if (CS%id_t_area > 0) then - tmp(:,:) = 1.0; call integrate_over_ice_sheet_area(G, ISS, tmp, 1.0, val) - call post_scalar_data(CS%id_t_area,val,CS%diag) - endif - if (CS%id_g_area > 0 .or. CS%id_f_area > 0) then - tmp(:,:) = 1.0; call masked_var_grounded(G,CS%dCS,tmp,tmp) - if (CS%id_g_area > 0) then - call integrate_over_ice_sheet_area(G, ISS, tmp, 1.0, val) - call post_scalar_data(CS%id_g_area,val,CS%diag) - endif - if (CS%id_f_area > 0) then - call integrate_over_ice_sheet_area(G, ISS, 1.0-tmp, 1.0, val) - call post_scalar_data(CS%id_f_area,val,CS%diag) - endif - endif + call process_and_post_scalar_data(CS, vaf0, vaf0_A, vaf0_G, Itime_step, dh_adott, dh_bdott) call disable_averaging(CS%diag) call cpu_clock_end(id_clock_shelf) @@ -926,20 +874,43 @@ subroutine shelf_calc_flux(sfc_state_in, fluxes_in, Time, time_step_in, CS) end subroutine shelf_calc_flux -subroutine integrate_over_ice_sheet_area(G, ISS, var, var_scale, var_out) +subroutine integrate_over_ice_sheet_area(G, ISS, var, var_scale, var_out, hemisphere) type(ocean_grid_type), intent(in) :: G !< The grid structure used by the ice shelf. type(ice_shelf_state), intent(in) :: ISS !< A structure with elements that describe the ice-shelf state real, dimension(SZI_(G),SZJ_(G)), intent(in) :: var !< Ice variable to integrate in arbitrary units [A ~> a] real, intent(in) :: var_scale !< Dimensional scaling for variable to integrate [a A-1 ~> 1] real, intent(out) :: var_out !< Variable integrated over the area of the ice sheet in arbitrary units [a m2] + integer, optional, intent(in) :: hemisphere !< 0 for Antarctica only, 1 for Greenland only. Otherwise, all ice sheets + integer :: IS_ID ! local copy of hemisphere real, dimension(SZI_(G),SZJ_(G)) :: var_cell !< Variable integrated over the ice-sheet area of each cell !! in arbitrary units [a m2] + integer, dimension(SZI_(G),SZJ_(G)) :: mask ! a mask for active cells depending on hemisphere indicated integer :: i,j + if (present(hemisphere)) then + IS_ID=hemisphere + else + IS_ID=-1 + endif + + mask(:,:)=0 + if (IS_ID==0) then !Antarctica (S. Hemisphere) only + do j = G%jsc,G%jec; do i = G%isc,G%iec + if (ISS%hmask(i,j)>0 .and. G%geoLatT(i,j)<=0.0) mask(i,j)=1 + enddo; enddo + elseif (IS_ID==1) then !Greenland (N. Hemisphere) only + do j = G%jsc,G%jec; do i = G%isc,G%iec + if (ISS%hmask(i,j)>0 .and. G%geoLatT(i,j)>0.0) mask(i,j)=1 + enddo; enddo + else !All ice sheets + mask(G%isc:G%iec,G%jsc:G%jec)=ISS%hmask(G%isc:G%iec,G%jsc:G%jec) + endif + var_cell(:,:)=0.0 do j = G%jsc,G%jec; do i = G%isc,G%iec - if (ISS%hmask(i,j)>0) var_cell(i,j) = (var(i,j) * var_scale) * (ISS%area_shelf_h(i,j) * G%US%L_to_m**2) + if (mask(i,j)>0) var_cell(i,j) = (var(i,j) * var_scale) * (ISS%area_shelf_h(i,j) * G%US%L_to_m**2) enddo; enddo + var_out = reproducing_sum(var_cell) end subroutine integrate_over_ice_sheet_area @@ -2031,11 +2002,12 @@ subroutine initialize_ice_shelf(param_file, ocn_grid, Time, CS, diag, Time_init, 'ice shelf surface mass flux deposition from atmosphere', & 'kg m-2 s-1', conversion=US%RZ_T_to_kg_m2s) endif - !scalars (area integrated) + + !scalars (area integrated over all ice sheets) CS%id_vaf = register_scalar_field('ice_shelf_model', 'int_vaf', CS%diag%axesT1, CS%Time, & 'Area integrated ice sheet volume above floatation', 'm3') CS%id_adott = register_scalar_field('ice_shelf_model', 'int_a', CS%diag%axesT1, CS%Time, & - 'Area integrated (entire ice sheet) change in ice-sheet thickness ' //& + 'Area integrated change in ice-sheet thickness ' //& 'due to surface accum+melt during a DT_THERM time step', 'm3') CS%id_g_adott = register_scalar_field('ice_shelf_model', 'int_a_ground', CS%diag%axesT1, CS%Time, & 'Area integrated change in grounded ice-sheet thickness ' //& @@ -2051,16 +2023,16 @@ subroutine initialize_ice_shelf(param_file, ocn_grid, Time, CS, diag, Time_init, CS%id_bdott_accum = register_scalar_field('ice_shelf_model', 'int_b_accum', CS%diag%axesT1, CS%Time, & 'Area integrated basal accumulation over ice shelves during a DT_THERM a time step', 'm3') CS%id_t_area = register_scalar_field('ice_shelf_model', 'tot_area', CS%diag%axesT1, CS%Time, & - 'Total area of entire ice-sheet', 'm2') + 'Total ice-sheet area', 'm2') CS%id_f_area = register_scalar_field('ice_shelf_model', 'tot_area_float', CS%diag%axesT1, CS%Time, & 'Total area of floating ice shelves', 'm2') CS%id_g_area = register_scalar_field('ice_shelf_model', 'tot_area_ground', CS%diag%axesT1, CS%Time, & - 'Total area of grounded ice sheet', 'm2') - !scalars (area integrated rates) + 'Total area of grounded ice sheets', 'm2') + !scalars (area integrated rates over all ice sheets) CS%id_dvafdt = register_scalar_field('ice_shelf_model', 'int_vafdot', CS%diag%axesT1, CS%Time, & 'Area integrated rate of change in ice-sheet volume above floatation', 'm3 s-1') CS%id_adot = register_scalar_field('ice_shelf_model', 'int_adot', CS%diag%axesT1, CS%Time, & - 'Area integrated (full ice sheet) rate of change in ice-sheet thickness due to surface accum+melt', 'm3 s-1') + 'Area integrated rate of change in ice-sheet thickness due to surface accum+melt', 'm3 s-1') CS%id_g_adot = register_scalar_field('ice_shelf_model', 'int_adot_ground', CS%diag%axesT1, CS%Time, & 'Area integrated rate of change in grounded ice-sheet thickness due to surface accum+melt', 'm3 s-1') CS%id_f_adot = register_scalar_field('ice_shelf_model', 'int_adot_float', CS%diag%axesT1, CS%Time, & @@ -2072,6 +2044,111 @@ subroutine initialize_ice_shelf(param_file, ocn_grid, Time, CS, diag, Time_init, CS%id_bdot_accum = register_scalar_field('ice_shelf_model', 'int_bdot_accum', CS%diag%axesT1, CS%Time, & 'Area integrated basal accumulation rate over ice shelves', 'm3 s-1') + !scalars (area integrated over the Antarctic ice sheet) + CS%id_Ant_vaf = register_scalar_field('ice_shelf_model', 'int_vaf_A', CS%diag%axesT1, CS%Time, & + 'Area integrated Antarctic ice sheet volume above floatation', 'm3') + CS%id_Ant_adott = register_scalar_field('ice_shelf_model', 'int_a_A', CS%diag%axesT1, CS%Time, & + 'Area integrated (Antarctic ice sheet) change in ice-sheet thickness ' //& + 'due to surface accum+melt during a DT_THERM time step', 'm3') + CS%id_Ant_g_adott = register_scalar_field('ice_shelf_model', 'int_a_ground_A', CS%diag%axesT1, CS%Time, & + 'Area integrated change in Antarctic grounded ice-sheet thickness ' //& + 'due to surface accum+melt during a DT_THERM time step', 'm3') + CS%id_Ant_f_adott = register_scalar_field('ice_shelf_model', 'int_a_float_A', CS%diag%axesT1, CS%Time, & + 'Area integrated change in Antarctic floating ice-shelf thickness ' //& + 'due to surface accum+melt during a DT_THERM time step', 'm3') + CS%id_Ant_bdott = register_scalar_field('ice_shelf_model', 'int_b_A', CS%diag%axesT1, CS%Time, & + 'Area integrated change in Antarctic floating ice-shelf thickness '//& + 'due to basal accum+melt during a DT_THERM time step', 'm3') + CS%id_Ant_bdott_melt = register_scalar_field('ice_shelf_model', 'int_b_melt_A', CS%diag%axesT1, CS%Time, & + 'Area integrated basal melt over Antarctic ice shelves during a DT_THERM time step', 'm3') + CS%id_Ant_bdott_accum = register_scalar_field('ice_shelf_model', 'int_b_accum_A', CS%diag%axesT1, CS%Time, & + 'Area integrated basal accumulation over Antarctic ice shelves during a DT_THERM a time step', 'm3') + CS%id_Ant_t_area = register_scalar_field('ice_shelf_model', 'tot_area_A', CS%diag%axesT1, CS%Time, & + 'Total area of Antarctic ice sheet', 'm2') + CS%id_Ant_f_area = register_scalar_field('ice_shelf_model', 'tot_area_float_A', CS%diag%axesT1, CS%Time, & + 'Total area of Antarctic floating ice shelves', 'm2') + CS%id_Ant_g_area = register_scalar_field('ice_shelf_model', 'tot_area_ground_A', CS%diag%axesT1, CS%Time, & + 'Total area of Antarctic grounded ice sheet', 'm2') + !scalars (area integrated rates over the Antarctic ice sheet) + CS%id_Ant_dvafdt = register_scalar_field('ice_shelf_model', 'int_vafdot_A', CS%diag%axesT1, CS%Time, & + 'Area integrated rate of change in Antarctic ice-sheet volume above floatation', 'm3 s-1') + CS%id_Ant_adot = register_scalar_field('ice_shelf_model', 'int_adot_A', CS%diag%axesT1, CS%Time, & + 'Area integrated rate of change in Antarctic ice-sheet thickness due to surface accum+melt', 'm3 s-1') + CS%id_Ant_g_adot = register_scalar_field('ice_shelf_model', 'int_adot_ground_A', CS%diag%axesT1, CS%Time, & + 'Area integrated rate of change in Antarctic grounded ice-sheet thickness due to surface accum+melt', 'm3 s-1') + CS%id_Ant_f_adot = register_scalar_field('ice_shelf_model', 'int_adot_float_A', CS%diag%axesT1, CS%Time, & + 'Area integrated rate of change in Antarctic floating ice-shelf thickness due to surface accum+melt', 'm3 s-1') + CS%id_Ant_bdot = register_scalar_field('ice_shelf_model', 'int_bdot_A', CS%diag%axesT1, CS%Time, & + 'Area integrated rate of change in Antarctic ice-shelf thickness due to basal accum+melt', 'm3 s-1') + CS%id_Ant_bdot_melt = register_scalar_field('ice_shelf_model', 'int_bdot_melt_A', CS%diag%axesT1, CS%Time, & + 'Area integrated basal melt rate over Antarctic ice shelves', 'm3 s-1') + CS%id_Ant_bdot_accum = register_scalar_field('ice_shelf_model', 'int_bdot_accum_A', CS%diag%axesT1, CS%Time, & + 'Area integrated basal accumulation rate over Antarctic ice shelves', 'm3 s-1') + + !scalars (area integrated over the Greenland ice sheet) + CS%id_Gr_vaf = register_scalar_field('ice_shelf_model', 'int_vaf_G', CS%diag%axesT1, CS%Time, & + 'Area integrated Greenland ice sheet volume above floatation', 'm3') + CS%id_Gr_adott = register_scalar_field('ice_shelf_model', 'int_a_G', CS%diag%axesT1, CS%Time, & + 'Area integrated (Greenland ice sheet) change in ice-sheet thickness ' //& + 'due to surface accum+melt during a DT_THERM time step', 'm3') + CS%id_Gr_g_adott = register_scalar_field('ice_shelf_model', 'int_a_ground_G', CS%diag%axesT1, CS%Time, & + 'Area integrated change in Greenland grounded ice-sheet thickness ' //& + 'due to surface accum+melt during a DT_THERM time step', 'm3') + CS%id_Gr_f_adott = register_scalar_field('ice_shelf_model', 'int_a_float_G', CS%diag%axesT1, CS%Time, & + 'Area integrated change in Greenland floating ice-shelf thickness ' //& + 'due to surface accum+melt during a DT_THERM time step', 'm3') + CS%id_Gr_bdott = register_scalar_field('ice_shelf_model', 'int_b_G', CS%diag%axesT1, CS%Time, & + 'Area integrated change in Greenland floating ice-shelf thickness '//& + 'due to basal accum+melt during a DT_THERM time step', 'm3') + CS%id_Gr_bdott_melt = register_scalar_field('ice_shelf_model', 'int_b_melt_G', CS%diag%axesT1, CS%Time, & + 'Area integrated basal melt over Greenland ice shelves during a DT_THERM time step', 'm3') + CS%id_Gr_bdott_accum = register_scalar_field('ice_shelf_model', 'int_b_accum_G', CS%diag%axesT1, CS%Time, & + 'Area integrated basal accumulation over Greenland ice shelves during a DT_THERM a time step', 'm3') + CS%id_Gr_t_area = register_scalar_field('ice_shelf_model', 'tot_area_G', CS%diag%axesT1, CS%Time, & + 'Total area of Greenland ice sheet', 'm2') + CS%id_Gr_f_area = register_scalar_field('ice_shelf_model', 'tot_area_float_G', CS%diag%axesT1, CS%Time, & + 'Total area of Greenland floating ice shelves', 'm2') + CS%id_Gr_g_area = register_scalar_field('ice_shelf_model', 'tot_area_ground_G', CS%diag%axesT1, CS%Time, & + 'Total area of Greenland grounded ice sheet', 'm2') + !scalars (area integrated rates over the Greenland ice sheet) + CS%id_Gr_dvafdt = register_scalar_field('ice_shelf_model', 'int_vafdot_G', CS%diag%axesT1, CS%Time, & + 'Area integrated rate of change in Greenland ice-sheet volume above floatation', 'm3 s-1') + CS%id_Gr_adot = register_scalar_field('ice_shelf_model', 'int_adot_G', CS%diag%axesT1, CS%Time, & + 'Area integrated rate of change in Greenland ice-sheet thickness due to surface accum+melt', 'm3 s-1') + CS%id_Gr_g_adot = register_scalar_field('ice_shelf_model', 'int_adot_ground_G', CS%diag%axesT1, CS%Time, & + 'Area integrated rate of change in Greenland grounded ice-sheet thickness due to surface accum+melt', 'm3 s-1') + CS%id_Gr_f_adot = register_scalar_field('ice_shelf_model', 'int_adot_float_G', CS%diag%axesT1, CS%Time, & + 'Area integrated rate of change in Greenland floating ice-shelf thickness due to surface accum+melt', 'm3 s-1') + CS%id_Gr_bdot = register_scalar_field('ice_shelf_model', 'int_bdot_G', CS%diag%axesT1, CS%Time, & + 'Area integrated rate of change in Greenland ice-shelf thickness due to basal accum+melt', 'm3 s-1') + CS%id_Gr_bdot_melt = register_scalar_field('ice_shelf_model', 'int_bdot_melt_G', CS%diag%axesT1, CS%Time, & + 'Area integrated basal melt rate over Greenland ice shelves', 'm3 s-1') + CS%id_Gr_bdot_accum = register_scalar_field('ice_shelf_model', 'int_bdot_accum_G', CS%diag%axesT1, CS%Time, & + 'Area integrated basal accumulation rate over Greenland ice shelves', 'm3 s-1') + + !Flags to calculate diagnostics related to surface/basal mass balance + if (CS%id_adott>0 .or. CS%id_g_adott>0 .or. CS%id_f_adott>0 .or. & + CS%id_adot >0 .or. CS%id_g_adot >0 .or. CS%id_f_adot >0 .or. & + CS%id_Ant_adott>0 .or. CS%id_Ant_g_adott>0 .or. CS%id_Ant_f_adott>0 .or. & + CS%id_Ant_adot >0 .or. CS%id_Ant_g_adot >0 .or. CS%id_Ant_f_adot >0 .or. & + CS%id_Gr_adott>0 .or. CS%id_Gr_g_adott>0 .or. CS%id_Gr_f_adott>0 .or. & + CS%id_Gr_adot >0 .or. CS%id_Gr_g_adot >0 .or. CS%id_Gr_f_adot >0) then + CS%smb_diag=.true. + else + CS%smb_diag=.false. + endif + + if (CS%id_bdott>0 .or. CS%id_bdott_melt>0 .or. CS%id_bdott_accum>0 .or. & + CS%id_bdot >0 .or. CS%id_bdot_melt >0 .or. CS%id_bdot_accum >0 .or. & + CS%id_Ant_bdott>0 .or. CS%id_Ant_bdott_melt>0 .or. CS%id_Ant_bdott_accum>0 .or. & + CS%id_Ant_bdot >0 .or. CS%id_Ant_bdot_melt >0 .or. CS%id_Ant_bdot_accum >0 .or. & + CS%id_Gr_bdott>0 .or. CS%id_Gr_bdott_melt>0 .or. CS%id_Gr_bdott_accum>0 .or. & + CS%id_Gr_bdot >0 .or. CS%id_Gr_bdot_melt >0 .or. CS%id_Gr_bdot_accum >0) then + CS%bmb_diag=.true. + else + CS%bmb_diag=.false. + endif + call MOM_IS_diag_mediator_close_registration(CS%diag) if (present(fluxes_in)) call initialize_ice_shelf_fluxes(CS, ocn_grid, US, fluxes_in) @@ -2447,11 +2524,9 @@ subroutine solo_step_ice_shelf(CS, time_interval, nsteps, Time, min_time_step_in logical :: coupled_GL ! If true the grounding line position is determined based on ! coupled ice-ocean dynamics. integer :: is, ie, js, je, i, j - real :: vaf0, vaf ! The previous and current volume above floatation [m3] - logical :: smb_diag=.false. ! Flags to calculate diagnostics related to surface/basal mass balance - real :: val ! Temporary value when calculating scalar diagnostics [various] + real :: vaf0, vaf0_A, vaf0_G !The previous volumes above floatation + !for all ice sheets, Antarctica only, or Greenland only [m3] real, dimension(SZI_(CS%grid),SZJ_(CS%grid)) :: & - tmp, & ! Temporary field used when calculating diagnostics [various] dh_adott_sum, & ! Surface melt/accumulation over a full time step, used for diagnostics [Z ~> m] dh_adott ! Surface melt/accumulation over a partial time step, used for diagnostics [Z ~> m] @@ -2475,14 +2550,14 @@ subroutine solo_step_ice_shelf(CS, time_interval, nsteps, Time, min_time_step_in ISS%dhdt_shelf(:,:) = ISS%h_shelf(:,:) - if (CS%id_adott>0 .or. CS%id_g_adott>0 .or. CS%id_f_adott>0 .or. & - CS%id_adot >0 .or. CS%id_g_adot >0 .or. CS%id_f_adot >0) then - smb_diag=.true. - dh_adott(:,:) = 0.0 ; dh_adott_sum(:,:) = 0.0 ; tmp(:,:) = 0.0 - endif + dh_adott(:,:)=0.0 + + if (CS%smb_diag) dh_adott_sum(:,:) = 0.0 - if (CS%id_dvafdt > 0) & !calculate previous volume above floatation - call volume_above_floatation(CS%dCS, G, ISS, vaf0) + !calculate previous volumes above floatation + if (CS%id_dvafdt > 0) call volume_above_floatation(CS%dCS, G, ISS, vaf0) !all ice sheet + if (CS%id_Ant_dvafdt > 0) call volume_above_floatation(CS%dCS, G, ISS, vaf0_A, hemisphere=0) !Antarctica only + if (CS%id_Gr_dvafdt > 0) call volume_above_floatation(CS%dCS, G, ISS, vaf0_G, hemisphere=1) !Greenland only do while (remaining_time > 0.0) nsteps = nsteps+1 @@ -2497,9 +2572,9 @@ subroutine solo_step_ice_shelf(CS, time_interval, nsteps, Time, min_time_step_in call MOM_mesg("solo_step_ice_shelf: "//mesg, 5) endif - if (smb_diag) dh_adott(is:ie,js:je) = ISS%h_shelf(is:ie,js:je) + if (CS%smb_diag) dh_adott(is:ie,js:je) = ISS%h_shelf(is:ie,js:je) call change_thickness_using_precip(CS, ISS, G, US, fluxes_in, time_step, Time) - if (smb_diag) dh_adott_sum(is:ie,js:je) = dh_adott_sum(is:ie,js:je) + & + if (CS%smb_diag) dh_adott_sum(is:ie,js:je) = dh_adott_sum(is:ie,js:je) + & (ISS%h_shelf(is:ie,js:je) - dh_adott(is:ie,js:je)) remaining_time = remaining_time - time_step @@ -2525,47 +2600,220 @@ subroutine solo_step_ice_shelf(CS, time_interval, nsteps, Time, min_time_step_in if (CS%id_h_shelf > 0) call post_data(CS%id_h_shelf ,ISS%h_shelf ,CS%diag) if (CS%id_dhdt_shelf > 0) call post_data(CS%id_dhdt_shelf ,ISS%dhdt_shelf ,CS%diag) if (CS%id_h_mask > 0) call post_data(CS%id_h_mask ,ISS%hmask ,CS%diag) - if (CS%id_vaf > 0 .or. CS%id_dvafdt > 0) & !calculate current volume above floatation (vaf) + call process_and_post_scalar_data(CS, vaf0, vaf0_A, vaf0_G, Ifull_time_step, dh_adott, dh_adott*0.0) + call disable_averaging(CS%diag) + + call IS_dynamics_post_data(full_time_step, Time, CS%dCS, G) +end subroutine solo_step_ice_shelf + +subroutine process_and_post_scalar_data(CS, vaf0, vaf0_A, vaf0_G, Itime_step, dh_adott, dh_bdott) + type(ice_shelf_CS), pointer :: CS !< A pointer to the ice shelf control structure + real :: vaf !< The current ice-sheet volume above floatation [m3] + real :: vaf0 !< The previous volumes above floatation for all ice sheets [m3] + real :: vaf0_A !< The previous volumes above floatation for the Antarctic ice sheet [m3] + real :: vaf0_G !< The previous volumes above floatation for the Greenland ice sheet [m3] + real :: Itime_step !< Inverse of the time step [T-1 ~> s-1] + real, dimension(SZI_(CS%grid),SZJ_(CS%grid)) :: & + dh_adott, & !< Surface (plus basal if solo shelf mode) melt/accumulation over a time step [Z ~> m] + dh_bdott, & !< Surface (plus basal if solo shelf mode) melt/accumulation over a time step [Z ~> m] + tmp ! Temporary field used when calculating diagnostics [various] + real :: val ! Temporary value when calculating scalar diagnostics [various] + type(ocean_grid_type), pointer :: G => NULL() ! A pointer to the ocean's grid structure + type(unit_scale_type), pointer :: US => NULL() ! Pointer to a structure containing various unit conversion factors + type(ice_shelf_state), pointer :: ISS => NULL() !< A structure with elements that describe the ice-shelf state + integer :: is, ie, js, je, i, j + + G => CS%grid + US => CS%US + ISS => CS%ISS + is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec + + !---ALL ICE SHEET---! + if (CS%id_vaf > 0 .or. CS%id_dvafdt > 0) & !calculate current volume above floatation (vaf) call volume_above_floatation(CS%dCS, G, ISS, vaf) - if (CS%id_vaf > 0) call post_scalar_data(CS%id_vaf ,vaf ,CS%diag) !current vaf - if (CS%id_dvafdt > 0) call post_scalar_data(CS%id_dvafdt,(vaf-vaf0)*Ifull_time_step,CS%diag) !d(vaf)/dt + if (CS%id_vaf > 0) call post_scalar_data(CS%id_vaf ,vaf ,CS%diag) !current vaf + if (CS%id_dvafdt > 0) call post_scalar_data(CS%id_dvafdt,(vaf-vaf0)*Itime_step,CS%diag) !d(vaf)/dt if (CS%id_adott > 0 .or. CS%id_adot > 0) then !surface accumulation - surface melt - call integrate_over_ice_sheet_area(G, ISS, dh_adott_sum, US%Z_to_m, val) - if (CS%id_adott > 0) call post_scalar_data(CS%id_adott,val ,CS%diag) - if (CS%id_adot > 0) call post_scalar_data(CS%id_adot ,val*Ifull_time_step,CS%diag) + call integrate_over_ice_sheet_area(G, ISS, dh_adott, US%Z_to_m, val) + if (CS%id_adott > 0) call post_scalar_data(CS%id_adott,val ,CS%diag) + if (CS%id_adot > 0) call post_scalar_data(CS%id_adot ,val*Itime_step,CS%diag) endif if (CS%id_g_adott > 0 .or. CS%id_g_adot > 0) then !grounded only: surface accumulation - surface melt - call masked_var_grounded(G,CS%dCS,dh_adott_sum,tmp) + call masked_var_grounded(G,CS%dCS,dh_adott,tmp) call integrate_over_ice_sheet_area(G, ISS, tmp, US%Z_to_m, val) - if (CS%id_g_adott > 0) call post_scalar_data(CS%id_g_adott,val ,CS%diag) - if (CS%id_g_adot > 0) call post_scalar_data(CS%id_g_adot ,val*Ifull_time_step,CS%diag) + if (CS%id_g_adott > 0) call post_scalar_data(CS%id_g_adott,val ,CS%diag) + if (CS%id_g_adot > 0) call post_scalar_data(CS%id_g_adot ,val*Itime_step,CS%diag) endif if (CS%id_f_adott > 0 .or. CS%id_f_adot > 0) then !floating only: surface accumulation - surface melt - call masked_var_grounded(G,CS%dCS,dh_adott_sum,tmp) - tmp(:,:) = dh_adott_sum(:,:) - tmp(:,:) + call masked_var_grounded(G,CS%dCS,dh_adott,tmp) + tmp(:,:) = dh_adott(:,:) - tmp(:,:) + call integrate_over_ice_sheet_area(G, ISS, tmp, US%Z_to_m, val) + if (CS%id_f_adott > 0) call post_scalar_data(CS%id_f_adott,val ,CS%diag) + if (CS%id_f_adot > 0) call post_scalar_data(CS%id_f_adot ,val*Itime_step,CS%diag) + endif + if (CS%id_bdott > 0 .or. CS%id_bdot > 0) then !bottom accumulation - bottom melt + call integrate_over_ice_sheet_area(G, ISS, dh_bdott, US%Z_to_m, val) + if (CS%id_bdott > 0) call post_scalar_data(CS%id_bdott,val ,CS%diag) + if (CS%id_bdot > 0) call post_scalar_data(CS%id_bdot ,val*Itime_step,CS%diag) + endif + if (CS%id_bdott_melt > 0 .or. CS%id_bdot_melt > 0) then !bottom melt + tmp(:,:)=0.0 + do j=js,je ; do i=is,ie + if (dh_bdott(i,j) < 0) tmp(i,j) = -dh_bdott(i,j) + enddo; enddo call integrate_over_ice_sheet_area(G, ISS, tmp, US%Z_to_m, val) - if (CS%id_f_adott > 0) call post_scalar_data(CS%id_f_adott,val ,CS%diag) - if (CS%id_f_adot > 0) call post_scalar_data(CS%id_f_adot ,val*Ifull_time_step,CS%diag) + if (CS%id_bdott_melt > 0) call post_scalar_data(CS%id_bdott_melt,val ,CS%diag) + if (CS%id_bdot_melt > 0) call post_scalar_data(CS%id_bdot_melt ,val*Itime_step,CS%diag) endif - if (CS%id_t_area > 0) then + if (CS%id_bdott_accum > 0 .or. CS%id_bdot_accum > 0) then !bottom accumulation + tmp(:,:)=0.0 + do j=js,je ; do i=is,ie + if (dh_bdott(i,j) > 0) tmp(i,j) = dh_bdott(i,j) + enddo; enddo + call integrate_over_ice_sheet_area(G, ISS, tmp, US%Z_to_m, val) + if (CS%id_bdott_accum > 0) call post_scalar_data(CS%id_bdott_accum,val ,CS%diag) + if (CS%id_bdot_accum > 0) call post_scalar_data(CS%id_bdot_accum ,val*Itime_step,CS%diag) + endif + if (CS%id_t_area > 0) then !ice sheet area tmp(:,:) = 1.0; call integrate_over_ice_sheet_area(G, ISS, tmp, 1.0, val) call post_scalar_data(CS%id_t_area,val,CS%diag) endif if (CS%id_g_area > 0 .or. CS%id_f_area > 0) then tmp(:,:) = 1.0; call masked_var_grounded(G,CS%dCS,tmp,tmp) - if (CS%id_g_area > 0) then + if (CS%id_g_area > 0) then !grounded only ice sheet area call integrate_over_ice_sheet_area(G, ISS, tmp, 1.0, val) call post_scalar_data(CS%id_g_area,val,CS%diag) endif - if (CS%id_f_area > 0) then + if (CS%id_f_area > 0) then !floating only ice sheet area (ice shelf area) call integrate_over_ice_sheet_area(G, ISS, 1.0-tmp, 1.0, val) call post_scalar_data(CS%id_f_area,val,CS%diag) endif endif - call disable_averaging(CS%diag) - call IS_dynamics_post_data(full_time_step, Time, CS%dCS, G) -end subroutine solo_step_ice_shelf + !---ANTARCTICA ONLY---! + if (CS%id_Ant_vaf > 0 .or. CS%id_Ant_dvafdt > 0) & !calculate current volume above floatation (vaf) + call volume_above_floatation(CS%dCS, G, ISS, vaf, hemisphere=0) + if (CS%id_Ant_vaf > 0) call post_scalar_data(CS%id_Ant_vaf ,vaf ,CS%diag) !current vaf + if (CS%id_Ant_dvafdt > 0) call post_scalar_data(CS%id_Ant_dvafdt,(vaf-vaf0_A)*Itime_step,CS%diag) !d(vaf)/dt + if (CS%id_Ant_adott > 0 .or. CS%id_Ant_adot > 0) then !surface accumulation - surface melt + call integrate_over_ice_sheet_area(G, ISS, dh_adott, US%Z_to_m, val, hemisphere=0) + if (CS%id_Ant_adott > 0) call post_scalar_data(CS%id_Ant_adott,val ,CS%diag) + if (CS%id_Ant_adot > 0) call post_scalar_data(CS%id_Ant_adot ,val*Itime_step,CS%diag) + endif + if (CS%id_Ant_g_adott > 0 .or. CS%id_Ant_g_adot > 0) then !grounded only: surface accumulation - surface melt + call masked_var_grounded(G,CS%dCS,dh_adott,tmp) + call integrate_over_ice_sheet_area(G, ISS, tmp, US%Z_to_m, val, hemisphere=0) + if (CS%id_Ant_g_adott > 0) call post_scalar_data(CS%id_Ant_g_adott,val ,CS%diag) + if (CS%id_Ant_g_adot > 0) call post_scalar_data(CS%id_Ant_g_adot ,val*Itime_step,CS%diag) + endif + if (CS%id_Ant_f_adott > 0 .or. CS%id_Ant_f_adot > 0) then !floating only: surface accumulation - surface melt + call masked_var_grounded(G,CS%dCS,dh_adott,tmp) + tmp(:,:) = dh_adott(:,:) - tmp(:,:) + call integrate_over_ice_sheet_area(G, ISS, tmp, US%Z_to_m, val, hemisphere=0) + if (CS%id_Ant_f_adott > 0) call post_scalar_data(CS%id_Ant_f_adott,val ,CS%diag) + if (CS%id_Ant_f_adot > 0) call post_scalar_data(CS%id_Ant_f_adot ,val*Itime_step,CS%diag) + endif + if (CS%id_Ant_bdott > 0 .or. CS%id_Ant_bdot > 0) then !bottom accumulation - bottom melt + call integrate_over_ice_sheet_area(G, ISS, dh_bdott, US%Z_to_m, val, hemisphere=0) + if (CS%id_Ant_bdott > 0) call post_scalar_data(CS%id_Ant_bdott,val ,CS%diag) + if (CS%id_Ant_bdot > 0) call post_scalar_data(CS%id_Ant_bdot ,val*Itime_step,CS%diag) + endif + if (CS%id_Ant_bdott_melt > 0 .or. CS%id_Ant_bdot_melt > 0) then !bottom melt + tmp(:,:)=0.0 + do j=js,je ; do i=is,ie + if (dh_bdott(i,j) < 0) tmp(i,j) = -dh_bdott(i,j) + enddo; enddo + call integrate_over_ice_sheet_area(G, ISS, tmp, US%Z_to_m, val, hemisphere=0) + if (CS%id_Ant_bdott_melt > 0) call post_scalar_data(CS%id_Ant_bdott_melt,val ,CS%diag) + if (CS%id_Ant_bdot_melt > 0) call post_scalar_data(CS%id_Ant_bdot_melt ,val*Itime_step,CS%diag) + endif + if (CS%id_Ant_bdott_accum > 0 .or. CS%id_Ant_bdot_accum > 0) then !bottom accumulation + tmp(:,:)=0.0 + do j=js,je ; do i=is,ie + if (dh_bdott(i,j) > 0) tmp(i,j) = dh_bdott(i,j) + enddo; enddo + call integrate_over_ice_sheet_area(G, ISS, tmp, US%Z_to_m, val, hemisphere=0) + if (CS%id_Ant_bdott_accum > 0) call post_scalar_data(CS%id_Ant_bdott_accum,val ,CS%diag) + if (CS%id_Ant_bdot_accum > 0) call post_scalar_data(CS%id_Ant_bdot_accum ,val*Itime_step,CS%diag) + endif + if (CS%id_Ant_t_area > 0) then !ice sheet area + tmp(:,:) = 1.0; call integrate_over_ice_sheet_area(G, ISS, tmp, 1.0, val, hemisphere=0) + call post_scalar_data(CS%id_Ant_t_area,val,CS%diag) + endif + if (CS%id_Ant_g_area > 0 .or. CS%id_Ant_f_area > 0) then + tmp(:,:) = 1.0; call masked_var_grounded(G,CS%dCS,tmp,tmp) + if (CS%id_Ant_g_area > 0) then !grounded only ice sheet area + call integrate_over_ice_sheet_area(G, ISS, tmp, 1.0, val, hemisphere=0) + call post_scalar_data(CS%id_Ant_g_area,val,CS%diag) + endif + if (CS%id_Ant_f_area > 0) then !floating only ice sheet area (ice shelf area) + call integrate_over_ice_sheet_area(G, ISS, 1.0-tmp, 1.0, val, hemisphere=0) + call post_scalar_data(CS%id_Ant_f_area,val,CS%diag) + endif + endif + + !---GREENLAND ONLY---! + if (CS%id_Gr_vaf > 0 .or. CS%id_Gr_dvafdt > 0) & !calculate current volume above floatation (vaf) + call volume_above_floatation(CS%dCS, G, ISS, vaf, hemisphere=1) + if (CS%id_Gr_vaf > 0) call post_scalar_data(CS%id_Gr_vaf ,vaf ,CS%diag) !current vaf + if (CS%id_Gr_dvafdt > 0) call post_scalar_data(CS%id_Gr_dvafdt,(vaf-vaf0_A)*Itime_step,CS%diag) !d(vaf)/dt + if (CS%id_Gr_adott > 0 .or. CS%id_Gr_adot > 0) then !surface accumulation - surface melt + call integrate_over_ice_sheet_area(G, ISS, dh_adott, US%Z_to_m, val, hemisphere=1) + if (CS%id_Gr_adott > 0) call post_scalar_data(CS%id_Gr_adott,val ,CS%diag) + if (CS%id_Gr_adot > 0) call post_scalar_data(CS%id_Gr_adot ,val*Itime_step,CS%diag) + endif + if (CS%id_Gr_g_adott > 0 .or. CS%id_Gr_g_adot > 0) then !grounded only: surface accumulation - surface melt + call masked_var_grounded(G,CS%dCS,dh_adott,tmp) + call integrate_over_ice_sheet_area(G, ISS, tmp, US%Z_to_m, val, hemisphere=1) + if (CS%id_Gr_g_adott > 0) call post_scalar_data(CS%id_Gr_g_adott,val ,CS%diag) + if (CS%id_Gr_g_adot > 0) call post_scalar_data(CS%id_Gr_g_adot ,val*Itime_step,CS%diag) + endif + if (CS%id_Gr_f_adott > 0 .or. CS%id_Gr_f_adot > 0) then !floating only: surface accumulation - surface melt + call masked_var_grounded(G,CS%dCS,dh_adott,tmp) + tmp(:,:) = dh_adott(:,:) - tmp(:,:) + call integrate_over_ice_sheet_area(G, ISS, tmp, US%Z_to_m, val, hemisphere=1) + if (CS%id_Gr_f_adott > 0) call post_scalar_data(CS%id_Gr_f_adott,val ,CS%diag) + if (CS%id_Gr_f_adot > 0) call post_scalar_data(CS%id_Gr_f_adot ,val*Itime_step,CS%diag) + endif + if (CS%id_Gr_bdott > 0 .or. CS%id_Gr_bdot > 0) then !bottom accumulation - bottom melt + call integrate_over_ice_sheet_area(G, ISS, dh_bdott, US%Z_to_m, val, hemisphere=1) + if (CS%id_Gr_bdott > 0) call post_scalar_data(CS%id_Gr_bdott,val ,CS%diag) + if (CS%id_Gr_bdot > 0) call post_scalar_data(CS%id_Gr_bdot ,val*Itime_step,CS%diag) + endif + if (CS%id_Gr_bdott_melt > 0 .or. CS%id_Gr_bdot_melt > 0) then !bottom melt + tmp(:,:)=0.0 + do j=js,je ; do i=is,ie + if (dh_bdott(i,j) < 0) tmp(i,j) = -dh_bdott(i,j) + enddo; enddo + call integrate_over_ice_sheet_area(G, ISS, tmp, US%Z_to_m, val, hemisphere=1) + if (CS%id_Gr_bdott_melt > 0) call post_scalar_data(CS%id_Gr_bdott_melt,val ,CS%diag) + if (CS%id_Gr_bdot_melt > 0) call post_scalar_data(CS%id_Gr_bdot_melt ,val*Itime_step,CS%diag) + endif + if (CS%id_Gr_bdott_accum > 0 .or. CS%id_Gr_bdot_accum > 0) then !bottom accumulation + tmp(:,:)=0.0 + do j=js,je ; do i=is,ie + if (dh_bdott(i,j) > 0) tmp(i,j) = dh_bdott(i,j) + enddo; enddo + call integrate_over_ice_sheet_area(G, ISS, tmp, US%Z_to_m, val, hemisphere=1) + if (CS%id_Gr_bdott_accum > 0) call post_scalar_data(CS%id_Gr_bdott_accum,val ,CS%diag) + if (CS%id_Gr_bdot_accum > 0) call post_scalar_data(CS%id_Gr_bdot_accum ,val*Itime_step,CS%diag) + endif + if (CS%id_Gr_t_area > 0) then !ice sheet area + tmp(:,:) = 1.0; call integrate_over_ice_sheet_area(G, ISS, tmp, 1.0, val, hemisphere=1) + call post_scalar_data(CS%id_Gr_t_area,val,CS%diag) + endif + if (CS%id_Gr_g_area > 0 .or. CS%id_Gr_f_area > 0) then + tmp(:,:) = 1.0; call masked_var_grounded(G,CS%dCS,tmp,tmp) + if (CS%id_Gr_g_area > 0) then !grounded only ice sheet area + call integrate_over_ice_sheet_area(G, ISS, tmp, 1.0, val, hemisphere=1) + call post_scalar_data(CS%id_Gr_g_area,val,CS%diag) + endif + if (CS%id_Gr_f_area > 0) then !floating only ice sheet area (ice shelf area) + call integrate_over_ice_sheet_area(G, ISS, 1.0-tmp, 1.0, val, hemisphere=1) + call post_scalar_data(CS%id_Gr_f_area,val,CS%diag) + endif + endif +end subroutine process_and_post_scalar_data !> \namespace mom_ice_shelf !! diff --git a/src/ice_shelf/MOM_ice_shelf_dynamics.F90 b/src/ice_shelf/MOM_ice_shelf_dynamics.F90 index c9a3baa8e2..15012353b1 100644 --- a/src/ice_shelf/MOM_ice_shelf_dynamics.F90 +++ b/src/ice_shelf/MOM_ice_shelf_dynamics.F90 @@ -955,26 +955,48 @@ subroutine update_ice_shelf(CS, ISS, G, US, time_step, Time, calve_ice_shelf_ber end subroutine update_ice_shelf -subroutine volume_above_floatation(CS, G, ISS, vaf) +subroutine volume_above_floatation(CS, G, ISS, vaf, hemisphere) type(ice_shelf_dyn_CS), intent(in) :: CS !< The ice shelf dynamics control structure type(ocean_grid_type), intent(in) :: G !< The grid structure used by the ice shelf. type(ice_shelf_state), intent(in) :: ISS !< A structure with elements that describe !! the ice-shelf state real, intent(out) :: vaf !< area integrated volume above floatation [m3] + integer, optional, intent(in) :: hemisphere !< 0 for Antarctica only, 1 for Greenland only. Otherwise, all ice sheets + integer :: IS_ID ! local copy of hemisphere real, dimension(SZI_(G),SZJ_(G)) :: vaf_cell !< cell-wise volume above floatation [m3] + integer, dimension(SZI_(G),SZJ_(G)) :: mask ! a mask for active cells depending on hemisphere indicated integer :: is,ie,js,je,i,j real :: rhoi_rhow, rhow_rhoi if (CS%GL_couple) & call MOM_error(FATAL, "MOM_ice_shelf_dyn, volume above floatation calculation assumes GL_couple=.FALSE..") - vaf_cell(:,:)=0.0 rhoi_rhow = CS%density_ice / CS%density_ocean_avg rhow_rhoi = CS%density_ocean_avg / CS%density_ice is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec + if (present(hemisphere)) then + IS_ID=hemisphere + else + IS_ID=-1 + endif + + mask(:,:)=0 + if (IS_ID==0) then !Antarctica (S. Hemisphere) only + do j = js,je; do i = is,ie + if (ISS%hmask(i,j)>0 .and. G%geoLatT(i,j)<=0.0) mask(i,j)=1 + enddo; enddo + elseif (IS_ID==1) then !Greenland (N. Hemisphere) only + do j = js,je; do i = is,ie + if (ISS%hmask(i,j)>0 .and. G%geoLatT(i,j)>0.0) mask(i,j)=1 + enddo; enddo + else !All ice sheets + mask(is:ie,js:je)=ISS%hmask(is:ie,js:je) + endif + + vaf_cell(:,:)=0.0 do j = js,je; do i = is,ie - if (ISS%hmask(i,j)>0) then + if (mask(i,j)>0) then if (CS%bed_elev(i,j) <= 0) then !grounded above sea level vaf_cell(i,j)= (ISS%h_shelf(i,j) * G%US%Z_to_m) * (ISS%area_shelf_h(i,j) * G%US%L_to_m**2) From d1322b26c06cfb08e1c2c3c881707a104eeb204c Mon Sep 17 00:00:00 2001 From: alex-huth Date: Mon, 19 Aug 2024 11:17:49 -0400 Subject: [PATCH 4/6] Added sx_shelf and sy_shelf as ice shelf diagnostics to save the surface slope fields used in the shallow shelf approximation (SSA). This is particularly helpful for determining whether unrealistic velocities are caused by unrealistically steep surface slopes, which can sometimes arise for example, on coarse grid cells that cover both a steep mountainous region and a realively flat ice shelf. Then, the MAX_SURFACE_SLOPE parameter can be tuned to set an upper bound on the SSA surface slope to avoid these steep-slope-induced problematic velocities. --- src/ice_shelf/MOM_ice_shelf_dynamics.F90 | 21 ++++++++++++++++++--- 1 file changed, 18 insertions(+), 3 deletions(-) diff --git a/src/ice_shelf/MOM_ice_shelf_dynamics.F90 b/src/ice_shelf/MOM_ice_shelf_dynamics.F90 index 15012353b1..d7c4de80b4 100644 --- a/src/ice_shelf/MOM_ice_shelf_dynamics.F90 +++ b/src/ice_shelf/MOM_ice_shelf_dynamics.F90 @@ -50,10 +50,14 @@ module MOM_ice_shelf_dynamics !! on q-points (B grid) [L T-1 ~> m s-1] real, pointer, dimension(:,:) :: v_shelf => NULL() !< the meridional velocity of the ice shelf/sheet !! on q-points (B grid) [L T-1 ~> m s-1] - real, pointer, dimension(:,:) :: taudx_shelf => NULL() !< the driving stress of the ice shelf/sheet + real, pointer, dimension(:,:) :: taudx_shelf => NULL() !< the zonal driving stress of the ice shelf/sheet !! on q-points (C grid) [R L2 T-2 ~> Pa] - real, pointer, dimension(:,:) :: taudy_shelf => NULL() !< the meridional stress of the ice shelf/sheet + real, pointer, dimension(:,:) :: taudy_shelf => NULL() !< the meridional driving stress of the ice shelf/sheet !! on q-points (C grid) [R L2 T-2 ~> Pa] + real, pointer, dimension(:,:) :: sx_shelf => NULL() !< the zonal surface slope of the ice shelf/sheet + !! on q-points (B grid) [nondim] + real, pointer, dimension(:,:) :: sy_shelf => NULL() !< the meridional surface slope of the ice shelf/sheet + !! on q-points (B grid) [nondim] real, pointer, dimension(:,:) :: u_face_mask => NULL() !< mask for velocity boundary conditions on the C-grid !! u-face - this is because the FEM cares about FACES THAT GET INTEGRATED OVER, !! not vertices. Will represent boundary conditions on computational boundary @@ -226,7 +230,8 @@ module MOM_ice_shelf_dynamics integer :: id_u_shelf = -1, id_v_shelf = -1, id_t_shelf = -1, & id_taudx_shelf = -1, id_taudy_shelf = -1, id_bed_elev = -1, & id_ground_frac = -1, id_col_thick = -1, id_OD_av = -1, id_float_cond = -1, & - id_u_mask = -1, id_v_mask = -1, id_ufb_mask =-1, id_vfb_mask = -1, id_t_mask = -1 + id_u_mask = -1, id_v_mask = -1, id_ufb_mask =-1, id_vfb_mask = -1, id_t_mask = -1, & + id_sx_shelf = -1, id_sy_shelf = -1 !>@} ! ids for outputting intermediate thickness in advection subroutine (debugging) !>@{ Diagnostic handles for debugging @@ -353,6 +358,8 @@ subroutine register_ice_shelf_dyn_restarts(G, US, param_file, CS, restart_CS) allocate(CS%ground_frac(isd:ied,jsd:jed), source=0.0) allocate(CS%taudx_shelf(IsdB:IedB,JsdB:JedB), source=0.0) allocate(CS%taudy_shelf(IsdB:IedB,JsdB:JedB), source=0.0) + allocate(CS%sx_shelf(isd:ied,jsd:jed), source=0.0) + allocate(CS%sy_shelf(isd:ied,jsd:jed), source=0.0) allocate(CS%bed_elev(isd:ied,jsd:jed), source=0.0) allocate(CS%u_bdry_val(IsdB:IedB,JsdB:JedB), source=0.0) allocate(CS%v_bdry_val(IsdB:IedB,JsdB:JedB), source=0.0) @@ -802,6 +809,10 @@ subroutine initialize_ice_shelf_dyn(param_file, Time, ISS, CS, G, US, diag, new_ 'x-driving stress of ice', 'kPa', conversion=1.e-3*US%RLZ_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-3*US%RLZ_T2_to_Pa) + CS%id_sx_shelf = register_diag_field('ice_shelf_model','sx_shelf',CS%diag%axesB1, Time, & + 'x-surface slope of ice', 'none') + CS%id_sy_shelf = register_diag_field('ice_shelf_model','sy_shelf',CS%diag%axesB1, Time, & + 'y-surface slope of ice', 'none') 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, & @@ -1047,6 +1058,8 @@ subroutine IS_dynamics_post_data(time_step, Time, CS, G) taud_y(:,:) = CS%taudy_shelf(:,:)*G%IareaBu(:,:) call post_data(CS%id_taudy_shelf, taud_y, CS%diag) endif + if (CS%id_sx_shelf > 0) call post_data(CS%id_sx_shelf, CS%sx_shelf, CS%diag) + if (CS%id_sy_shelf > 0) call post_data(CS%id_sy_shelf, CS%sy_shelf, CS%diag) if (CS%id_ground_frac > 0) call post_data(CS%id_ground_frac, CS%ground_frac, CS%diag) if (CS%id_float_cond > 0) call post_data(CS%id_float_cond, CS%float_cond, CS%diag) if (CS%id_OD_av >0) call post_data(CS%id_OD_av, CS%OD_av,CS%diag) @@ -2441,6 +2454,8 @@ subroutine calc_shelf_driving_stress(CS, ISS, G, US, taudx, taudy, OD) sx_e(i,j) = (-.25 * G%areaT(i,j)) * ((rho * grav) * (max(ISS%h_shelf(i,j),CS%min_h_shelf) * sx)) sy_e(i,j) = (-.25 * G%areaT(i,j)) * ((rho * grav) * (max(ISS%h_shelf(i,j),CS%min_h_shelf) * sy)) + CS%sx_shelf(i,j) = sx ; CS%sy_shelf(i,j) = sy + !Stress (Neumann) boundary conditions if (CS%ground_frac(i,j) == 1) then neumann_val = ((.5 * grav) * (rho * max(ISS%h_shelf(i,j),CS%min_h_shelf)**2 - rhow * CS%bed_elev(i,j)**2)) From 419f1cd2ec3e5d76f9073b0af7b9bdc76b318569 Mon Sep 17 00:00:00 2001 From: alex-huth Date: Mon, 19 Aug 2024 12:05:43 -0400 Subject: [PATCH 5/6] Fix doxygen errors for ice-sheet process_and_post_scalar_data routine --- src/ice_shelf/MOM_ice_shelf.F90 | 16 +++++++++------- 1 file changed, 9 insertions(+), 7 deletions(-) diff --git a/src/ice_shelf/MOM_ice_shelf.F90 b/src/ice_shelf/MOM_ice_shelf.F90 index 18e3669f98..52b1ebdcea 100644 --- a/src/ice_shelf/MOM_ice_shelf.F90 +++ b/src/ice_shelf/MOM_ice_shelf.F90 @@ -2606,21 +2606,23 @@ subroutine solo_step_ice_shelf(CS, time_interval, nsteps, Time, min_time_step_in call IS_dynamics_post_data(full_time_step, Time, CS%dCS, G) end subroutine solo_step_ice_shelf +!> Post_data calls for ice-sheet scalars subroutine process_and_post_scalar_data(CS, vaf0, vaf0_A, vaf0_G, Itime_step, dh_adott, dh_bdott) type(ice_shelf_CS), pointer :: CS !< A pointer to the ice shelf control structure - real :: vaf !< The current ice-sheet volume above floatation [m3] real :: vaf0 !< The previous volumes above floatation for all ice sheets [m3] real :: vaf0_A !< The previous volumes above floatation for the Antarctic ice sheet [m3] real :: vaf0_G !< The previous volumes above floatation for the Greenland ice sheet [m3] real :: Itime_step !< Inverse of the time step [T-1 ~> s-1] - real, dimension(SZI_(CS%grid),SZJ_(CS%grid)) :: & - dh_adott, & !< Surface (plus basal if solo shelf mode) melt/accumulation over a time step [Z ~> m] - dh_bdott, & !< Surface (plus basal if solo shelf mode) melt/accumulation over a time step [Z ~> m] - tmp ! Temporary field used when calculating diagnostics [various] - real :: val ! Temporary value when calculating scalar diagnostics [various] + real, dimension(SZI_(CS%grid),SZJ_(CS%grid)) :: dh_adott !< Surface (plus basal if solo shelf mode) + !! melt/accumulation over a time step [Z ~> m] + real, dimension(SZI_(CS%grid),SZJ_(CS%grid)) :: dh_bdott !< Surface (plus basal if solo shelf mode) + !! melt/accumulation over a time step [Z ~> m] + real, dimension(SZI_(CS%grid),SZJ_(CS%grid)) :: tmp ! Temporary field used when calculating diagnostics [various] + real :: vaf ! The current ice-sheet volume above floatation [m3] + real :: val ! Temporary value when calculating scalar diagnostics [various] type(ocean_grid_type), pointer :: G => NULL() ! A pointer to the ocean's grid structure type(unit_scale_type), pointer :: US => NULL() ! Pointer to a structure containing various unit conversion factors - type(ice_shelf_state), pointer :: ISS => NULL() !< A structure with elements that describe the ice-shelf state + type(ice_shelf_state), pointer :: ISS => NULL() ! A structure with elements that describe the ice-shelf state integer :: is, ie, js, je, i, j G => CS%grid From 7c18900c351283f06b06e6396fd187d9330861b6 Mon Sep 17 00:00:00 2001 From: alex-huth Date: Tue, 17 Sep 2024 15:49:50 -0400 Subject: [PATCH 6/6] FMA fix associated with enforcement of max allowed ice-shelf surface slope --- src/ice_shelf/MOM_ice_shelf_dynamics.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/ice_shelf/MOM_ice_shelf_dynamics.F90 b/src/ice_shelf/MOM_ice_shelf_dynamics.F90 index d7c4de80b4..3ed262e5f3 100644 --- a/src/ice_shelf/MOM_ice_shelf_dynamics.F90 +++ b/src/ice_shelf/MOM_ice_shelf_dynamics.F90 @@ -2447,7 +2447,7 @@ subroutine calc_shelf_driving_stress(CS, ISS, G, US, taudx, taudy, OD) endif if (CS%max_surface_slope>0) then - scale = min(CS%max_surface_slope/sqrt(sx**2+sy**2),1.0) + scale = min(CS%max_surface_slope/sqrt((sx**2)+(sy**2)),1.0) sx = scale*sx; sy = scale*sy endif