OCN to GLC thermal forcing coupling

cime_config/tests.py

@@ -266,6 +266,7 @@
             "SMS_D_Ld1.T62_oQU240wLI.GMPAS-IAF-PISMF.mpaso-impl_top_drag",
             "SMS_D_Ld1.T62_oQU240.GMPAS-IAF.mpaso-harmonic_mean_drag",
             "SMS_D_Ld1.T62_oQU240.GMPAS-IAF.mpaso-upwind_advection",
+            "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

@@ -80,6 +80,8 @@ module glc_comp_mct
 
   integer :: glcLogUnit ! unit number for glc log
 
+  logical  :: ocn_c2_glctf           ! .true.  => ocn to glc thermal forcing coupling on
+
   ! MPAS Datatypes
   !type (dm_info), pointer :: dminfo
   type (core_type), pointer :: corelist => null()
@@ -522,6 +524,9 @@ end subroutine xml_stream_get_attributes
     ! Determine coupling type (not currently needed by MALI)
     call seq_infodata_GetData(infodata, cpl_seq_option=cpl_seq_option)
 
+    ! Determine if ocn to glc thermal forcing coupling is on
+    call seq_infodata_GetData(infodata, ocn_c2_glctf=ocn_c2_glctf)
+
     ! Initialize the MALI core
     ierr = domain % core % core_init(domain, timeStamp)
     if ( ierr /= 0 ) then
@@ -1383,7 +1388,8 @@ subroutine glc_import_mct(x2g_g, errorCode)
                                                floatingBasalMassBal,&
                                                surfaceTemperature,&
                                                basalOceanHeatflx,&
-                                               OceanDensity
+                                               OceanDensity, &
+                                               ismip6_2dThermalForcing
 
     errorCode = 0
 
@@ -1401,13 +1407,16 @@ 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)
+            if (ocn_c2_glctf) &
+               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

@@ -354,6 +354,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',perrWith='quiet')
 
     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

@@ -108,6 +108,8 @@ module ocn_comp_mct
 
    integer  :: nsend, nrecv
 
+   logical  :: ocn_c2_glctf           ! .true.  => ocn to glc thermal forcing coupling on
+
    character(len=StrKIND) :: runtype, coupleTimeStamp
 
    type(seq_infodata_type), pointer :: infodata
@@ -223,6 +225,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
@@ -304,6 +307,9 @@ end subroutine xml_stream_get_attributes
     ! Determine coupling type
     call seq_infodata_GetData(infodata, cpl_seq_option=cpl_seq_option)
 
+    ! Determine if ocn to glc thermal forcing coupling is on
+    call seq_infodata_GetData(infodata, ocn_c2_glctf=ocn_c2_glctf)
+
 !-----------------------------------------------------------------------
 !
 !   initialize the model run
@@ -878,6 +884,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 +2705,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 +2725,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 +2736,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 +2792,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 +2956,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' .and. ocn_c2_glctf) 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."