Merge branch 'matthewhoffman/mali/gis_ocn_to_glc_coupling' into next …

…(PR #6632)

OCN to GLC thermal forcing coupling

This pull request introduces coupling from OCN to GLC to pass 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. Facemelting as a function of
ocean thermal state is a critical OCN/GLC coupling for Greenland, but it
is relevant to Antarctica as well. This OCN-GLC thermal forcing coupling
is activated whenever OCN is present and GLC is active.

This PR has pieces in MPAS-Ocean, the coupler, and MALI:
* In MPAS-Ocean, a new avgThermalForcingAtCritDepth field is added and
  calculated. The depth at which to calculate the thermal forcing is
  namelist configurable. The default is 300 m.
* In MALI, the field passed from the coupler is collected in
  ismip6_2dThermalForcing and used in the existing facemelting
  parameterization. Namelist and streams are updated to use facemelting
  in all simulations and output relevant variables. For compsets where
  TF is not being coupled, it will have values of 0 and result in no
  facemelting being applied; it is safe to have this option generally
  turned on.
* In the coupler, a new remapping object is created that handles the TF
  mapping independently of any other fields passed from OCN to GLC. New
  mapping files are added. Any coupling variables related to the
  existing ice-shelf coupling now include the 'shelf' suffix and all
  coupling variables related to this new thermal-forcing coupling have
  the 'tf' suffix.
* A new compset is created that is a standard G-case but with active
  MALI. New mapping files are added between the relevant grids for the
  new mapper.
* The new OCN2GLC_TF_SMAPNAME mapping should be a nearest source to
  destination mapping with the MPAS-Ocean source mesh masked to only
  include grid cells that are deeper than the critical depth to avoid
  extrapolating undefined values. A tool to generate that map file is
  MPAS-Tools
* Two new model_grids are added that are compatible with this feature
  and include the special mapping file required:
  TL319_IcoswISC30E3r5_gis1to10kmR2 and TL319_oQU240wLI_gis20
* The new coupling is controlled by a new namelist option in MPAS-Ocean:
  config_glc_thermal_forcing_coupling_mode. It defaults to 'off' and
  there are no compsets that explicitly set it for now, making this a
  stealth feature.

[NML] for new mappings, new ocean variable
[BFB] stealth feature (except A-cases due to new coupler field)

note: testing will report this as a FAIL due to:
"BASELINE master: FIELDLIST field lists differ (otherwise bit-for-bit)"

cime_config/tests.py

@@ -269,6 +269,7 @@
             "SMS_D_Ld1.T62_oQU240.GMPAS-IAF.mpaso-harmonic_mean_drag",
             "SMS_D_Ld1.T62_oQU240.GMPAS-IAF.mpaso-upwind_advection",
             "ERS_Ld5_D.T62_oQU240.GMPAS-IAF.mpaso-conservation_check",
+            "ERS_Ld5.TL319_oQU240wLI_gis20.MPAS_LISIO_JRA1p5.mpaso-ocn_glc_tf_coupling",
             )
         },
 

components/mpas-albany-landice/bld/namelist_files/namelist_defaults_mali.xml

@@ -98,7 +98,7 @@
 <config_temperature_profile_variability_period>1.0</config_temperature_profile_variability_period>
 <config_temperature_profile_variability_phase>0.0</config_temperature_profile_variability_phase>
 <config_temperature_profile_GL_depth_fraction>0.25</config_temperature_profile_GL_depth_fraction>
-<config_front_mass_bal_grounded>'none'</config_front_mass_bal_grounded>
+<config_front_mass_bal_grounded>'ismip6'</config_front_mass_bal_grounded>
 <config_use_3d_thermal_forcing_for_face_melt>.false.</config_use_3d_thermal_forcing_for_face_melt>
 <config_beta_ocean_thermal_forcing>1.18</config_beta_ocean_thermal_forcing>
 <config_add_ocean_thermal_forcing>0.0</config_add_ocean_thermal_forcing>

components/mpas-albany-landice/cime_config/buildnml

@@ -88,7 +88,7 @@ def buildnml(case, caseroot, compname):
         decomp_date += '051920'
         decomp_prefix += 'mpasli.graph.info.'
     elif glc_grid == 'mpas.gis1to10kmR2':
-        grid_date += '20230202'
+        grid_date += '20240513'
         grid_prefix += 'gis_1to10km_r02'
         decomp_date += '020223'
         decomp_prefix += 'mpasli.graph.info.'
@@ -247,6 +247,9 @@ def buildnml(case, caseroot, compname):
         lines.append('    <var name="calvingThickness"/>')
         lines.append('    <var name="restoreThickness"/>')
         lines.append('    <var name="dHdt"/>')
+        lines.append('    <var name="ismip6_2dThermalForcing"/>')
+        lines.append('    <var name="faceMeltSpeed"/>')
+        lines.append('    <var name="faceMeltingThickness"/>')
         lines.append('    <var name="deltat"/>')
         lines.append('    <var name="daysSinceStart"/>')
         lines.append('    <var name="simulationStartTime"/>')

components/mpas-albany-landice/driver/glc_comp_mct.F

@@ -1387,7 +1387,8 @@ subroutine glc_import_mct(x2g_g, errorCode)
                                                floatingBasalMassBal,&
                                                surfaceTemperature,&
                                                basalOceanHeatflx,&
-                                               OceanDensity
+                                               OceanDensity, &
+                                               ismip6_2dThermalForcing
 
     errorCode = 0
 
@@ -1405,13 +1406,15 @@ subroutine glc_import_mct(x2g_g, errorCode)
          call mpas_pool_get_array(geometryPool, 'sfcMassBal', sfcMassBal)
          call mpas_pool_get_array(geometryPool, 'floatingBasalMassBal',floatingBasalMassBal)
          call mpas_pool_get_array(thermalPool, 'surfaceTemperature',surfaceTemperature)
+         call mpas_pool_get_array(geometryPool, 'ismip6_2dThermalForcing', ismip6_2dThermalForcing)
 !         call mpas_pool_get_array(thermalPool, 'basalOceanHeatflx',basalOceanHeatflx)
          !call mpas_pool_get_array(geometryPool, 'OceanDensity',OceanDensity)
 
          do i = 1, nCellsSolve
             n = n + 1
             sfcMassBal(i) = x2g_g % rAttr(index_x2g_Flgl_qice, n)
             floatingBasalMassBal(i) = x2g_g % rAttr(index_x2g_Fogx_qiceli, n)
+            ismip6_2dThermalForcing(i) = x2g_g % rAttr(index_x2g_So_tf2d, n)
 !            surfaceTemperature(i) = x2g_g % rAttr(index_x2g_Sl_tsrf, n)
 !JW            basalOceanHeatflx(i) = x2g_g % rAttr(index_x2g_Fogo_qiceh, n)
 !            basalOceanHeatflx(i) = x2g_g % rAttr(index_x2g_Fogx_qicehi, n)

components/mpas-albany-landice/driver/glc_cpl_indices.F

@@ -22,6 +22,7 @@ module glc_cpl_indices
   integer, public :: index_x2g_So_htv    = 0 !Ice shelf ocean heat transfer velocity
   integer, public :: index_x2g_So_stv    = 0 !Ice shelf ocean salinity transfer velocity
   integer, public :: index_x2g_So_rhoeff = 0 !Ocean effective pressure
+  integer, public :: index_x2g_So_tf2d   = 0 !Ocean thermal forcing at predefined critical depth
   integer, public :: index_x2g_Fogx_qiceli = 0 !Subshelf mass flux
   integer, public :: index_x2g_Fogx_qicehi = 0 !Subshelf heat flux for the ice sheet
 
@@ -70,6 +71,7 @@ subroutine glc_cpl_indices_set( )
     index_x2g_Fogx_qiceli =               mct_avect_indexra(x2g,'Fogx_qiceli',perrwith='quiet')
     index_x2g_Fogx_qicehi =               mct_avect_indexra(x2g,'Fogx_qicehi',perrwith='quiet')
     index_x2g_So_rhoeff =                 mct_avect_indexra(x2g,'So_rhoeff',perrwith='quiet')
+    index_x2g_So_tf2d =                   mct_avect_indexra(x2g,'So_tf2d',perrwith='quiet')
 
     !Following block of x2g/g2x vectors are used internally within coupler for subshelf melt flux
     !calculations (and so do not have directly-related export-side arrays)

components/mpas-ocean/bld/build-namelist

@@ -731,6 +731,8 @@ if (($OCN_ICEBERG eq 'true') && ($OCN_FORCING eq 'active_atm')) {
 } else {
 	add_default($nl, 'config_remove_ais_ice_runoff', 'val'=>".false.");
 }
+add_default($nl, 'config_glc_thermal_forcing_coupling_mode');
+add_default($nl, 'config_2d_thermal_forcing_depth');
 
 ######################################
 # Namelist group: shortwaveRadiation #

components/mpas-ocean/bld/build-namelist-section

@@ -239,6 +239,8 @@ add_default($nl, 'config_sgr_salinity_prescribed');
 
 add_default($nl, 'config_remove_ais_river_runoff');
 add_default($nl, 'config_remove_ais_ice_runoff');
+add_default($nl, 'config_glc_thermal_forcing_coupling_mode');
+add_default($nl, 'config_2d_thermal_forcing_depth');
 
 ######################################
 # Namelist group: shortwaveRadiation #

components/mpas-ocean/bld/namelist_files/namelist_defaults_mpaso.xml

@@ -364,6 +364,8 @@
 <!-- coupling -->
 <config_remove_ais_river_runoff>.false.</config_remove_ais_river_runoff>
 <config_remove_ais_ice_runoff>.false.</config_remove_ais_ice_runoff>
+<config_glc_thermal_forcing_coupling_mode>'off'</config_glc_thermal_forcing_coupling_mode>
+<config_2d_thermal_forcing_depth>300.0</config_2d_thermal_forcing_depth>
 
 <!-- shortwaveRadiation -->
 <config_sw_absorption_type>'jerlov'</config_sw_absorption_type>

components/mpas-ocean/bld/namelist_files/namelist_definition_mpaso.xml

@@ -1271,6 +1271,22 @@ Valid values: .true. or .false.
 Default: Defined in namelist_defaults.xml
 </entry>
 
+<entry id="config_glc_thermal_forcing_coupling_mode" type="char*1024"
+	category="coupling" group="coupling">
+If and how MPAS-Ocean sends thermal forcing to GLC (MALI) in E3SM.  This is used for ocean coupling with a melt parameterization for grounded marine ice-cliffs in MALI.  This is primarily relevant to the Greenland Ice Sheet, but also relevant to the Antarctic Ice Sheet. 'none' means no coupling of thermal forcing.  '2d' means thermal forcing at a prescribed depth is passed to GLC.  That depth is controlled by 'config_2d_thermal_forcing_depth', and the resulting thermal forcing field is calculated in the field 'avgThermalForcingAtCritDepth'.
+
+Valid values: 'off', '2d'
+Default: Defined in namelist_defaults.xml
+</entry>
+
+<entry id="config_2d_thermal_forcing_depth" type="real"
+	category="coupling" group="coupling">
+Depth at which to pass 2d thermal forcing to the coupler for use in the GLC component.  Note that mapping files for this field must be created with a mask to exclude ocean grid cells shallower than this value and thus must be regenerated if this value is changed.
+
+Valid values: any non-negative value
+Default: Defined in namelist_defaults.xml
+</entry>
+
 
 <!-- shortwaveRadiation -->
 

components/mpas-ocean/cime_config/testdefs/testmods_dirs/mpaso/ocn_glc_tf_coupling/README

@@ -0,0 +1,12 @@
+This testdef is used to test a stealth feature that enables coupling between
+OCN and GLC for Greenland, which passes ocean thermal forcing from OCN to GLC
+and uses that in a parameterization for marine melting of grounded vertical
+cliffs.
+
+It changes one mpaso namelist variable,
+    config_glc_thermal_forcing_coupling_mode
+from its default value to '2d'.
+This tests the ocn/glc TF coupling.
+
+It also specified that DATM forcing should be restricted to 1958.
+This allows JRA1p5 forcing to be used without a large input data requirement.

components/mpas-ocean/cime_config/testdefs/testmods_dirs/mpaso/ocn_glc_tf_coupling/shell_commands

@@ -0,0 +1,4 @@
+./xmlchange DATM_CLMNCEP_YR_START=1958
+./xmlchange DATM_CLMNCEP_YR_END=1958
+./xmlchange DROF_STRM_YR_START=1958
+./xmlchange DROF_STRM_YR_END=1958

components/mpas-ocean/cime_config/testdefs/testmods_dirs/mpaso/ocn_glc_tf_coupling/user_nl_mpaso

@@ -0,0 +1 @@
+config_glc_thermal_forcing_coupling_mode = '2d'

components/mpas-ocean/driver/mpaso_cpl_indices.F

@@ -37,6 +37,7 @@ module mpaso_cpl_indices
   integer :: index_o2x_So_htv !ocean heat-transfer velocity
   integer :: index_o2x_So_stv !ocean salt-transfer velocity
   integer :: index_o2x_So_rhoeff !ocean effective density
+  integer :: index_o2x_So_tf2d !ocean thermal forcing at predefined critical depth
 
 
   ! ocn -> drv  (BGC)
@@ -208,6 +209,7 @@ subroutine mpaso_cpl_indices_set( )
     index_o2x_So_htv        = mct_avect_indexra(o2x,'So_htv')
     index_o2x_So_stv        = mct_avect_indexra(o2x,'So_stv')
     index_o2x_So_rhoeff     = mct_avect_indexra(o2x,'So_rhoeff')
+    index_o2x_So_tf2d       = mct_avect_indexra(o2x,'So_tf2d')
 
     index_o2x_So_algae1     = mct_avect_indexra(o2x,'So_algae1',perrWith='quiet')
     index_o2x_So_algae2     = mct_avect_indexra(o2x,'So_algae2',perrWith='quiet')

components/mpas-ocean/driver/ocn_comp_mct.F

@@ -223,6 +223,7 @@ subroutine ocn_init_mct( EClock, cdata_o, x2o_o, o2x_o, NLFilename )!{{{
     logical, pointer :: config_use_activeTracers_surface_restoring
     logical, pointer :: config_use_surface_salinity_monthly_restoring
     character (len=StrKIND), pointer :: config_land_ice_flux_mode
+    character (len=StrKIND), pointer :: config_glc_thermal_forcing_coupling_mode
 
     ! ssh coupling interval initialization
     integer, pointer :: index_avgZonalSSHGradient, index_avgMeridionalSSHGradient
@@ -878,6 +879,16 @@ end subroutine xml_stream_get_attributes
        call mpas_log_write('ERROR: unknown land_ice_flux_mode: ' // trim(config_land_ice_flux_mode), MPAS_LOG_CRIT)
     end if
 
+    call mpas_pool_get_config(domain % configs, 'config_glc_thermal_forcing_coupling_mode', config_glc_thermal_forcing_coupling_mode)
+    if ( trim(config_glc_thermal_forcing_coupling_mode) == 'off' ) then
+       call seq_infodata_PutData(infodata, ocn_c2_glctf=.false.)
+    else if ( trim(config_glc_thermal_forcing_coupling_mode) == '2d' ) then
+       call seq_infodata_PutData(infodata, ocn_c2_glctf=.true.)
+    else
+       call mpas_log_write('ERROR: unknown config_glc_thermal_forcing_coupling_mode: ' // &
+               trim(config_glc_thermal_forcing_coupling_mode), MPAS_LOG_CRIT)
+    end if
+
 !-----------------------------------------------------------------------
 !
 !   get initial state from driver
@@ -2689,7 +2700,8 @@ subroutine ocn_export_mct(o2x_o, errorCode)   !{{{
                                                avgRemovedRiverRunoffFlux, &
                                                avgRemovedIceRunoffFlux, &
                                                avgLandIceHeatFlux, &
-                                               avgRemovedIceRunoffHeatFlux
+                                               avgRemovedIceRunoffHeatFlux, &
+                                               avgThermalForcingAtCritDepth
 
    real (kind=RKIND), dimension(:,:), pointer :: avgTracersSurfaceValue, avgSurfaceVelocity, &
                                                  avgSSHGradient, avgOceanSurfacePhytoC, &
@@ -2708,6 +2720,7 @@ subroutine ocn_export_mct(o2x_o, errorCode)   !{{{
                        config_use_MacroMoleculesTracers_sea_ice_coupling
 
    character (len=StrKIND), pointer :: config_land_ice_flux_mode
+   character (len=StrKIND), pointer :: config_glc_thermal_forcing_coupling_mode
 
    logical :: keepFrazil
 
@@ -2718,6 +2731,8 @@ subroutine ocn_export_mct(o2x_o, errorCode)   !{{{
    call mpas_pool_get_config(domain % configs, 'config_land_ice_flux_mode', config_land_ice_flux_mode)
    call mpas_pool_get_config(domain % configs, 'config_remove_ais_river_runoff', config_remove_ais_river_runoff)
    call mpas_pool_get_config(domain % configs, 'config_remove_ais_ice_runoff', config_remove_ais_ice_runoff)
+   call mpas_pool_get_config(domain % configs, 'config_glc_thermal_forcing_coupling_mode', &
+                                                config_glc_thermal_forcing_coupling_mode)
    call mpas_pool_get_config(domain % configs, 'config_use_DMSTracers', config_use_DMSTracers)
    call mpas_pool_get_config(domain % configs, 'config_use_MacroMoleculesTracers', config_use_MacroMoleculesTracers)
    call mpas_pool_get_config(domain % configs, 'config_use_ecosysTracers_sea_ice_coupling',  &
@@ -2772,6 +2787,9 @@ subroutine ocn_export_mct(o2x_o, errorCode)   !{{{
         call mpas_pool_get_array(forcingPool, 'avgRemovedIceRunoffFlux', avgRemovedIceRunoffFlux)
         call mpas_pool_get_array(forcingPool, 'avgRemovedIceRunoffHeatFlux', avgRemovedIceRunoffHeatFlux)
      endif
+     if (trim(config_glc_thermal_forcing_coupling_mode) == '2d') then
+        call mpas_pool_get_array(forcingPool, 'avgThermalForcingAtCritDepth', avgThermalForcingAtCritDepth)
+     endif
 
      ! BGC fields
      if (config_use_ecosysTracers) then
@@ -2933,6 +2951,10 @@ subroutine ocn_export_mct(o2x_o, errorCode)   !{{{
           o2x_o % rAttr(index_o2x_So_stv, n) = landIceTracerTransferVelocities(indexSaltTrans,i)
           o2x_o % rAttr(index_o2x_So_rhoeff, n) = 0.0_RKIND
        endif
+       if (trim(config_glc_thermal_forcing_coupling_mode) == '2d') then
+          o2x_o % rAttr(index_o2x_So_tf2d, n) = avgThermalForcingAtCritDepth(i)
+       endif
+
 
        !Fyke: test
        !write(stderrUnit,*) 'n=',n

components/mpas-ocean/src/Registry.xml

@@ -807,6 +807,14 @@
 					description="If true, solid runoff from the Antarctic Ice Sheet (below 60S latitude) coming from the coupled is zeroed in the coupler import routines.  To be used with data iceberg fluxes coming from the sea ice model."
 					possible_values=".true. or .false."
 		/>
+		<nml_option name="config_glc_thermal_forcing_coupling_mode" type="character" default_value="off"
+					description="If and how MPAS-Ocean sends thermal forcing to GLC (MALI) in E3SM.  This is used for ocean coupling with a melt parameterization for grounded marine ice-cliffs in MALI.  This is primarily relevant to the Greenland Ice Sheet, but also relevant to the Antarctic Ice Sheet. 'none' means no coupling of thermal forcing.  '2d' means thermal forcing at a prescribed depth is passed to GLC.  That depth is controlled by 'config_2d_thermal_forcing_depth', and the resulting thermal forcing field is calculated in the field 'avgThermalForcingAtCritDepth'."
+					possible_values="'off', '2d'"
+		/>
+		<nml_option name="config_2d_thermal_forcing_depth" type="real" default_value="300" units="m"
+					description="Depth at which to pass 2d thermal forcing to the coupler for use in the GLC component.  Note that mapping files for this field must be created with a mask to exclude ocean grid cells shallower than this value and thus must be regenerated if this value is changed."
+					possible_values="any non-negative value"
+		/>
 	</nml_record>
 	<nml_record name="shortwaveRadiation" mode="init;forward">
 		<nml_option name="config_sw_absorption_type" type="character" default_value="none"
@@ -4010,6 +4018,9 @@
 			description="The time-averaged effective ocean density within ice shelves based on Archimedes' principle."
 			packages="landIceCouplingPKG"
 		/>
+		<var name="avgThermalForcingAtCritDepth" type="real" dimensions="nCells Time" units="C"
+			description="The time-averaged thermal forcing at the predefined critical depth specified by config_2d_thermal_forcing_depth"
+		/>
 		<!-- Input fields for data (prescribed) land-ice fluxes -->
 		<var name="dataLandIceFreshwaterFlux" type="real" dimensions="nCells Time" units="kg m^-2 s^-1"
 			description="Flux of mass through the ocean surface, as read in from a forcing file. Positive into ocean."