Ice Sheet Coupling

Overview

Ice sheet coupling with the rest of the model depends on which ice sheet is being modeled, Antarctica or Greenland.

Ice-shelf coupling: (ocn_c2_glcshelf, glcshelf_c2_ocn) Fluxes between the bottom of the ice shelf and the ocean underneath it. The ice-shelf coupling will be used for Antarctica (left side of figure).  Greenland has some small ice shelves, but they are too small to be resolved sufficiently in MPAS-Ocean, so we do not intend to use the ice-shelf coupling for Greenland cases. glcshelf coupling is unique in doing a lot of calculations inside the coupler.  It’s calculating the melt and heat fluxes on the ocean coupling interval (finer than GLC) and the GLC mesh (finer than OCN). 

Thermal-forcing coupling: (ocn_c2_glctf) Passes thermal forcing at a prescribed ocean depth (300 m) from MPAS-Ocean to MALI, where MALI uses it to calculate grounded marine melting through an existing 'facemelting' parameterization. The thermal-forcing coupling is primarily intended for Greenland where there are primarily vertical cliffs and little floating ice (right side of figure). In that geometry, MPAS-Ocean is not equipped to calculate melt rates, and even if it could, most of the narrow fjords around Greenland would be unresolved so the ocean domain rarely intercepts the Greenland marine termini.  Instead, we make use of marine melting parameterizations that are a function of the thermal forcing (ocean temperature minus in situ freezing temperature).  We also make use of a capability we added to MALI to horizontally extrapolate ocean thermal forcing from wherever on the MALI mesh it’s available to the current terminus positions.  While the thermal-forcing coupling is primarily developed for Greenland, Antarctica also has vertical marine cliffs without ice shelves in some places (and will likely have more in the future), so we do anticipate eventually using the TF coupling in Antarctica at the same time as the ice-shelf coupling. Implemented in https://github.com/E3SM-Project/E3SM/pull/6632

land coupling: (glc_c2_lnd and lnd_c2_glc) the surface mass balance (SMB) between the ice sheet and the atmosphere is calculated in the land model. The “glc_2_lnd” coupling passes the necessary data to and from that calculation. The land ice passes elevation of the ice sheet in a flexible number of topography classes (10 in the example below) to the SMB routine in the land model.

sea-ice coupling: (glc_c2_ice): icebergs. Icebergs form from the ice sheets but are then advected by the sea-ice model.

ocean coupling: (glc_c2_ocn): Sg_icemask, Sg_blit, Sg_blis, Sg_lithop, Sg_icemask_grounded

There is no ocn_c2_glc, ice_c2_glcshelf., ice_c2_glc

Attribute strings

land to glc adds the following coupling fields

• l2x states: Sl_tsrf00 through Sl_tsrf10  and Sl_topo00 through Sl_topo10

• l2x states to glc: Sl_tsrf00 to Sl_tsrf10 AND Sl_topo00 to Sl_topo10

  • Sl_tsrf: Surface temperature of glacier (one for each elevation class)

  • Sl_topo: surface height of glacier (one for each elevation class)

• l2x fluxes: Flgl_qice00 through Flgl_qice10

• l2x_fluxes_to_glc:  Flgl_qice00 to Flgl_qice10

  • Flgl_qice: New glacier ice flux

glc to lnd adds the following coupling fields:

• x2l states: Sg_ice_covered00 through Sg_ice_covered10  AND  Sg_topo00 through Sg_topo10

• x2l states from glc: Sg_ice_covered00 through Sg_ice_covered10 AND Sg_topo00 through Sg_topo10

  • Sg_ice_covered: Fraction of glacier area

  • Sg_topo: surface height of glacier

• g2x states to lnd: Sg_icemask:Sg_icemask_coupled_fluxes:Sg_ice_covered:Sg_topo

  • Sg_icemask: Ice sheet grid coverage on global grid

  • Sg_icemask_coupled_fluxes: Ice sheet mask where we are potentially sending non-zero fluxes

  • Sg_ice_covered: Fraction of glacier area

  • Sg_topo

• glc2lnd_non_ec_fields: Sg_icemask, Sg_icemask_coupled_fluxes

• x2l fluxes:  Flgg_hflx00 to Flgg_hflx10

• x2l fluxes from glc:  Flgg_hflx00 to Flgg_hflx10

  • Flgg_hflx: Downward heat flux from glacier interior.

• glc2lnd_ec_extra_fields: Flgg_hflx

glc to ocn coupling involves these fields

• g2x fluxes: Fogg_rofl, Fogg_rofi (both set to 0), Fogx_qicelo, Fogx_qiceho

• g2o_liq_fluxes: Fogg_rofl

• g2o_iceq_fluxes: Fogg_rofi

• x2o fluxes: Fogx_qicelo, Fogx_qiceho (CALCULATED IN COUPLER)

• x2o states: Sg_icemask, Sg_blit, Sg_blis, Sg_lithop, Sg_icemask_grounded, Sg_icemask_floating, Sg_tbot, Sg_dztbot

• x2g states: So_blt:So_bls:So_htv:So_stv:So_rhoeff (NONE USED)

• x2g fluxes: Fogx_qiceli:Fogx_qicehi

Resulting land ice states and fluxes

• g2x_states: Sg_icemask, Sg_icemask_coupled_fluxes, Sg_ice_covered, Sg_topo, Sg_blit, Sg_blis, Sg_lithop, Sg_icemask_grounded, Sg_icemask_floating, Sg_tbot, Sg_dztbot

  • Sg_blit: Boundary layer interface temperature for ocean (NOT SENT FROM MALI)

  • Sg_blis: Boundary layer interface salinity for ocean (NOT SENT FROM MALI)

  • Sg_lithop: Ice sheet lithostatic pressure

  • Sg_icemask_grounded: Grounded ice mask

  • Sg_icemask_floating: Floating ice mask

  • Sg_tbot: Bottom layer ice temperature

  • Sg_dztbot: Bottom layer ice layer half thickness

• g2x_fluxes: Fogg_rofl, Fogg_rofi, Figg_rofi, Flgg_hflx, Fogx_qicelo, Fogx_qiceho

  • Fogg_rofl glc liquid runoff flux to ocean (Set to 0 in MALI)

  • Fogg_rofi glc frozen runoff flux to ocean (Set to 0 in MALI)

  • Figg_rofi: glc frozen runoff_iceberg flux to sea ice (Set to 0 in MALI)

  • Flgg_hflx: Downward heat flux from glacier interior (NOT SENT FROM MALI)

  • Flgg_qicelo: Subshelf liquid flux for ocean (NOT SENT FROM MALI)

  • Fogx_qiceho: Subshelf heat flux for the ocean (NOT SENT FROM MALI)

  • g2o_liq_fluxes: Fogg_rofl

  • g2o_iceq_fluxes: Fogg_rofi

• x2g_states: Sl_tsrf, So_blt, So_bls, So_htv, So_stv, So_rhoeff (NONE USED IN MALI)

  • Sl_tsrf:

  • So_blt: Ice shelf boundary layer ocean temperature

  • So_bls: Ice shelf boundary layer ocean salinity

  • So_htv: Ice shelf ocean heat transfer velocity

  • So_stv: Ice shelf ocean salinity transfer velocity

  • So_rhoeff: Ocean effective pressure

  • x2g_states_from_lnd: Sl_tsrf

• x2g_fluxes: Flgl_qice, Fogx_qiceli, Fogx_qicehi

  • Flgl_qice: New glacier ice flux (USED IN MALI)

  • Fogx_qiceli: Subshelf mass flux for ice sheet (USED IN MALI)

  • Fogx_qicehi: Subshelf heat flux for ice sheet (NOT USED)

  • x2g_fluxes_from_lnd: Flgl_qice

Modifications to coupler flow.

Assuming glc_present and glc_prognostic are TRUE so not including.

init

prep_lnd_init (glc_c2_lnd)

if(glc_c2_lnd) init mapper_Sg2l, mapper_Fg2l; derive field lists glc2lnd_non_ec_fields and glc2lnd_ec_extra_fields

prep_ocn_init (glc_c2_ocn, glcshelf_c2_ocn)

if (glc_c2_ocn) init mapper_Rg2o_liq and mapper_Rg2o_ice

if(glcshelf_c2_ocn) init mapper_Sg2o and mapper_Fg2o

prep_ice_init (glc_c2_ice, glcshelf_c2_ice)

if (glc_c2_ice) init mapper_Rg2i

if(glcshelf_c2_ice) init mapper_Sg2i,mapper_Fg2i

prep_glc_init( lnd_c2_glc, ocn_c2_glctf, ocn_c2_glcshelf)

if ( lnd_c2_glc .or. do_hist_l2x1yrg) initialize l2gacc_lx with seq_flds_l2x_fields_to_glc

if ( lnd_c2_glc) initialize l2x_gx with seq_flds_x2g_fields

if (lnd_c2_glc) init mapper_Sl2g, mapper_Fl2g, mapper_Fg2l; form g2x_lx_fields

if (ocn_c2_glctf OR ocn_c2_glcshelf) init o2x_gx with seq_flds_o2x_fields; init x2gacc_gx

also init mapper_So2g, mapper_Fo2g and arrays for compute_melt_fluxes

if (ocn_c2_glctf) init mapper_So2g_tf

if(ocn_c2_glcshelf) init mapper_So2g_shelf

if (glc_c2_ocn) do mapping Rg2o_liq, Rg2o_ice; g2x_gx to g2x_ox; seq_flds_g2o_liq/ice_fluxes,

if (glcshelf_c2_ocn) do mapping Sg2o, Fg2o; g2x_gx to g2x_ox MIGHT OVERWRITE ABOVE

if(glc_c2_ice) do mapping Rg2i; g2x_gx to g2x_ix ‘Fixx_rofi’ only

if(glcshelf_c2_ice) do mapping_Sg2i; g2x_gx, g2x_ix, Sg_icemask_coupled_fluxes

if (glc_c2_lnd) prep_lnd_calc_g2x_lx do mapping with Fg2l, from g2x_gx to g2x_lx, only glc2lnd_non_ec_fields (Sg_icemask, Sg_icemask_coupled_fluxes)

call map_glc2lnd_ec do special mapping from glc to lnd elevation classes using mapper_Fg2l, glc2lnd_ec_extra_fields

in run loop

after land receive

call cime_run_glc_setup_send

if (lnd_c2_glc OR ocn_c2_glctf OR ocn_c2_glcshelf)

if (glcrun_avg_alarm) accumulate for lnd and ocn

if (lnd_c2_glc) call prep_glc_calc_l2x_gx

call prep_glc_map_qice_conservative_lnd2glc: special (bilinear with correction) mapping of Flgl_qice only. See Bilinear lnd to glc.

map_lnd2glc (mapper_Sl2g)

if (smb_renormalize) prep_glc_renormalize_smb;(calls map_glc2lnd_ec)

call prep_glc_map_one_state_field_lnd2glc: special bilinear mapping one state at a time: also call map_lnd2glc

if (ocn_c2_glctf) call prep_glc_mrg_ocn

else prep_glc_zero_flds

At bottom after running atm and ocn

run glc (pass in glc_prognostic)

call cime_run_glc_recv_post

component_exch(glc, flow='c2x')

call cime_run_ocn_recv_post

call cime_run_ocnglc_coupling()

if (ocn_c2_glctf OR (glcshelf_to_ocn AND ocn_to_glcshelf))

call prep_glc_calc_o2x_gx:

if (ocn_c2_gtctf) call mapper_So2g_tf; o2x_ox to o2x_gx seq_flds_x2g_tf_states_from_ocn

if (ocn_c2_glcshelf) call mapper_So2g_shelf; o2x_ox to o2x_gx seq_flds_x2g_shelf_states_from_ocn

if (ocn_c2_glcshelf AND glcshelf_c2_ocn)

call prep_glc_calculate_subshelf_boundary_fluxes

call mapper_So2g_shelf; o2x_ox, x2g_gx; seq_flds_x2g_shelf_states_from_ocn

compute_melt_fluxes

call prep_ocn_shelf_calc_g2x_ox: map_Sg2o, map_Fg2o; g2x_gx, g2x_ox

call prep_glc_accum_ocn (accumulate x2g_g fields)

if ( glcshelf_c2_ice) call mapper map_Sg2i; g2x_gx, g2x_ix; Sg_icemask_coupled_fluxes

Special functions in land ice coupling

prep_glc_accum_avg

called by cime_run_glc_setup_send

         call prep_glc_accum_avg(timer='CPL:glcprep_avg', &
              lnd2glc_averaged_now=lnd2glc_averaged_now)

accumulation Avs: l2gacc_lx and x2gacc_gx are local save private Avs in prep_glc_mod

if (l2gacc_lx_cnt > 1) call mct_avect_avg(l2gacc_lx(eli), l2gacc_lx_cnt)

if (x2gacc_gx_cnt > 1) then

call mct_avect_avg(x2gacc_gx(egi), x2gacc_gx_cnt)

call mct_avect_copy(x2gacc_gx(egi), x2g_gx)

prep_glc_calc_l2x_gx

called by cime_run_glc_setup_send

call prep_glc_calc_l2x_gx(fractions_lx, timer='CPL:glcprep_lnd2glc')

loop over 1, number of fields in(seq_flds_x2g_fluxes_from_lnd) which is 1: Flgl_qice

call prep_glc_map_qice_conservative_lnd2glc

end loop

loop over 1, number of fields in seq_flds_x2g_states_from_lnd which is 1: Sl_tsrf

call prep_glc_map_one_state_field_lnd2glc

end loop

prep_glc_map_qice_conservate_lnd2glc

called by prep_glc_calc_l2x_gx which is called by cime_run_glc_setup_send

    ! Use a bilinear (Sl2g) mapper, as for states.
    ! The Fg2l mapper is needed to map some glc fields to the land grid
    !  for purposes of conservation.
    call prep_glc_map_qice_conservative_lnd2glc(egi=egi, eli=eli, &
              fractions_lx = fractions_lx(efi), &
              mapper_Sl2g = mapper_Sl2g, &
              mapper_Fg2l = mapper_Fg2l)

using the accumulated l2gacc_lx, map only Flgl_qice (the SMB) using map_lnd2glc with mapper_Sl2g and l2x_gx

export Flgl_qice from l2x_gx to an array qice_g

multiply qice_g by area_g/aream_g

if (smb_renormalize) call prep_glc_renormalize_smb using g2x_gx with mapper_Fg2l. Also modifies qice_g

import qice_g back in to l2x_gx in the Flgl_qice attribute.

prep_glc_map_one_state_field_lnd2glc

called only by prep_glc_calc_l2x_gx which is called by cime_run_glc_setup_send

called only for Sl_tsrf

calls map_lnd2glc for Sl_tsrf mapping local l2gacc_lx to l2x_gx. Also uses fields from g2x_gx

map_glc2lnd_ec

called by prep_lnd_calc_g2x_lx and prep_glc_renormalize_smb

   call map_glc2lnd_ec( &
        g2x_g = g2x_gx, &
        frac_field = glc_frac_field, &
        topo_field = glc_topo_field, &
        icemask_field = glc_icemask_field, &
        extra_fields = glc2lnd_ec_extra_fields, &
        mapper = mapper_Fg2l, &
        g2x_l = g2x_lx(egi)

glc_frac_field = ‘Sg_ice_covered’

glc_topo_field = ‘Sg_topo’

glc_icemask_field = 'Sg_icemask'

glc2lnd_ec_extra_fields = Flgg_hflx

This will map Sg_ice_covered, Sg_topo and Flgg_hflx from glc mesh to land model mesh with elevation classes: Sg_ice_covered00, Sg_ice_covered01, etc. Sg_topo00, Sg_topo01, etc., Flgg_hflx00, Flgg_hflx01, etc. There is a loop over elevations classes.

The call from prep_glc_renormalize_smb does not include extra_fields.

map_lnd2glc

called by both prep_glc_map_qice_conservative_lnd2glc and prep_glc_map_one_state_field_lnd2glc

maps one field from land with elevation classes to “glacier”

Note: “fieldname” for the “conservative” is Flql_qice which is calculated on multiple elevation classes in the land model.

export Sg_ice_covered and Sg_topo into arrays.

call get_glc_elevation_classes(glc_ice_covered, glc_topo, glc_elevclass) which returns the integer array glc_elevclass that has the elevation class index that glc_topo is contained in. glc_ice_covered is either 0 (bare land) or 1 (fully ice).

call map_bare_land (l2x_l, landfrac_l, “fieldname”, mapper, array) creates an glacier av with fieldname_bare_land = fieldname // elevclass_as_string as the only attribute and maps the one field (which is on elev class 0) to that glacier av from l2x_l. Then pulls out the field into an array and passes it back.

That field is copied to data_g

call map_ice_covered(l2x_l, landfrac_l, “fieldname”, glc_topo, mapper, array) see below. returns an array callled data_g_ice_covered

Form the final data_g field with

where (glc_elevclass /= 0)
   data_g = data_g_ice_covered
end where

map_ice_covered

called only by map_lnd2glc

create a list of of land model attributes with elevation classes except for class 0. Result will be something like 'Flgl_qice01:Flgl_qice02: ... :Flgl_qice10:Sl_topo01:Sl_topo02: ... :Sltopo10'

Create a glacer-sized aV with above attributes.

map l2x_l to that aV with normalization and weights lfrac

loop over elevation classes and export the Flgl_qice and Sl_topo to a 2D array for each elevation class and point.

Create the field on the glacier mesh by looking at the elevation classes.

do i= 1, size_of_glacier_mesh

if topo(i) lt mapped_topo(i,1) then

the actual topography is lower then the height of the lowest elevation class so assign the field value of the lowest elevation class: data(i) = mapped_data(i, 1)

if topo(i) gt mapped_topo(i,nEC) then

the actual topography is higher then highest mapped elevation class so assign the field value of the highest elevation class: data(i) = mapped_data(i,nEC)

else need to do an interpolation with the 2 elevation classes nearest the actual topo.

do j = 2, nEC

if (topo(i) < mapped_topo(i,j) then

elev_lower = mapped_topo(i,j-1)

elev_upper = mapped_topo(i,j)

ediff = elev_upper - elev_lower

data(i) = mapped_data(i,j-1) * (elev-upper - topo(i))/ediff + mapped_data(i,j) * (topo(i)- elev_lower)/ediff

enddo

endif

enddo

return the data array.

prep_glc_renormalize_smb

called only by prep_glc_map_qice_conservate_lnd2glc

! Renormalizes surface mass balance (smb, here named qice_g) so that the global
! integral on the glc grid is equal to the global integral on the land grid.

qice_g is passed in as an array.

map Sg_icemask to land grid with normalization. (BUG note from Bill Sacks: I think we actually want norm = .false. here, but this requires some more thought. See https://github.com/ESMCI/cime/issues/1516)

Then pull out mapped field in to an array.

map Sg_ice_covered and Sg_topo to land grid and elevation classes with map_glc2lnd_ec. No extra fields as in other call to map_glc2lnd_ec.

loop over elevation classes

export the mapped Sg_ice_coverd and local time averaged Flgl_qice to 2D arrays sized by (land grid, elevation classes)

Find the local sum of negative and positive Flgl_qice with a double loop over land points and elevation classes.

loop over land points

effective_area = min(lfrac(n),Sg_icemask_l(n)) * aream_l(n) where lfrac is from fractions_lx

loop over elevation classes

if (qice >= 0) local_accum_on_land_grid = local_accum_on_land_grid + effective_area * frac_l(n,ec) * qice_l(n,ec) frac_l is the mapped Sg_ice_covered.

else local_ablat_on_land_grid = local_ablat_on_land_grid + effective_area * frac_l(n,ec) * qice_l(n,ec

end double loop

Use shr_mpi_sum to find the global sum global_accum_on_land_grid and global_ablat_on_land_grid from the local sums computed above NOT BFB with PROC COUNT. Then bcast global sum back to all procs

Loop over glacier points and find the local weighted sum of qice_g (the input array)

if(qice >=0) local_accum_on_glc_grid = local_accum_on_glc_grid + Sg_icemask_g(n) * aream_g(n) * qice_g(n)

else local_ablat_on_glc_grid = local_ablat_on_glc_grid + Sg_icemask_g(n) * aream_g(n) * qice_g(n)

end loop

use shr_mpi_sum again to find global sums global_accum_on_glc_grid, global_ablat_on_glc_grid NOT BFB with PROC COUNT. Then bcast global sum back to all procs.

Find the renormalization factor

if (global_ablat_on_glc_grid < 0.0_r8) then  ! negative by definition
   ablat_renorm_factor = global_ablat_on_land_grid / global_ablat_on_glc_grid

if (global_accum_on_glc_grid > 0.0_r8) then
    accum_renorm_factor = global_accum_on_land_grid / global_accum_on_glc_grid

Apply renormalization

do n = 1, lsize_g
   if (qice_g(n) >= 0.0_r8) then
      qice_g(n) = qice_g(n) * accum_renorm_factor
   else
      qice_g(n) = qice_g(n) * ablat_renorm_factor
   endif
enddo

return.

prep_glc_calculate_subshelf_boundary_fluxes

called by cime_run_ocnglc_coupling

On the ice sheet grid, calculate shelf boundary fluxes

Pull out input arrays from attribute vectors all on glacier grid

So_blt, So_bls, So_htv, So_stv. Sg_tbot, Sg_dztbot, Sg_lithop, Sg_icemask_floating

call compute_melt_fluxes

take output from above and place in attributes: Sg_blis, Sg_blit, Fogx_qiceho, Fogx_gicelo, Fogx_qicehi, Fogx_qiceli

compute_melt_fluxes

a physics calculation that works only with Fortran arrays (no Avs). called only by prep_glc_calculate_subshelf_boundary_fluxes

(Text in routine from Xylar) This routine computes melt fluxes (melt rate, temperature fluxes into the ice and the ocean, and salt flux) as well as the interface temperature and salinity. This routine expects an ice temperature in the bottom layer of ice and ocean temperature and salinity in the top ocean layer as well as the pressure at the ice/ocean interface.
The ocean heat and salt transfer velocities are determined based on observations of turbulent mixing rates in the under-ice boundary layer. They should be the product of the friction velocity and a (possibly spatially variable) non-dimensional transfer coefficient.
The iceTemperatureDistance (Sg_dztbot) is the distance between the location where the iceTemperature is supplied and the ice-ocean interface, used to compute a temperature gradient. The ice thermal conductivity, SHR_CONST_KAPPA_LAND_ICE, is zero for the freezing solution from Holland and Jenkins (1999) in which the ice is purely insulating.