Add MCT ice-shelf melt and interface-temperature coupling between MPAS-Ocean and MALI

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

@@ -914,9 +914,9 @@ subroutine glc_run_mct( EClock, cdata_g, x2g_g, g2x_g)!{{{
          ! Finally, set whence to latest_before so we have piecewise-constant forcing.
          ! Could add, e.g., linear interpolation later.
 
-         ! NOTE: If any input ("forcing") files here contain fields that the coupler also passes, 
-         ! the input file versions here will clobber the fields passed from the coupler. 
-         
+         ! NOTE: If any input ("forcing") files here contain fields that the coupler also passes,
+         ! the input file versions here will clobber the fields passed from the coupler.
+
          call mpas_stream_mgr_read(domain % streamManager, whence=MPAS_STREAM_LATEST_BEFORE, saveActualWhen = .true., ierr=err_tmp)
          err = ior(err, err_tmp)
          call mpas_stream_mgr_reset_alarms(domain % streamManager, direction=MPAS_STREAM_INPUT, ierr=err_tmp)
@@ -1410,18 +1410,24 @@ subroutine glc_import_mct(x2g_g, errorCode)
    real (kind=RKIND), dimension(:), pointer :: sfcMassBal,&
                                                floatingBasalMassBal,&
                                                surfaceTemperature,&
-                                               basalOceanHeatflx,&
-                                               OceanDensity
+                                               basalOceanHeatFlux,&
+                                               iceOceanInterfaceTemperature
+
    real (kind=RKIND), dimension(:,:), pointer :: ismip6shelfMelt_3dThermalForcing
    integer, dimension(:,:), pointer :: orig3dOceanMask
-   real(kind=RKIND), pointer :: config_invalid_value_TF
+   real(kind=RKIND), pointer :: config_invalid_value_TF, config_min_floating_frac, &
+                                config_invalid_value_ocn_intr_temp
    character(len=StrKIND), pointer :: config_smb_source
    real(kind=RKIND) :: fractionalMaskVal
 
     errorCode = 0
 
     call mpas_pool_get_config(domain % configs, 'config_invalid_value_TF', config_invalid_value_TF)
     call mpas_pool_get_config(domain % configs, 'config_smb_source', config_smb_source)
+    call mpas_pool_get_config(domain % configs, 'config_min_floating_frac', config_min_floating_frac)
+    call mpas_pool_get_config(domain % configs, 'config_min_floating_frac', config_min_floating_frac)
+    call mpas_pool_get_config(domain % configs, 'config_invalid_value_ocn_intr_temp', &
+                              config_invalid_value_ocn_intr_temp)
 
     n = 0
     block => domain % blocklist
@@ -1440,10 +1446,10 @@ 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(thermalPool, 'basalOceanHeatFlux', basalOceanHeatFlux)
          call mpas_pool_get_array(geometryPool, 'ismip6shelfMelt_3dThermalForcing', ismip6shelfMelt_3dThermalForcing)
          call mpas_pool_get_array(extrapOceanDataPool, 'orig3dOceanMask', orig3dOceanMask)
-!         call mpas_pool_get_array(thermalPool, 'basalOceanHeatflx',basalOceanHeatflx)
-         !call mpas_pool_get_array(geometryPool, 'OceanDensity',OceanDensity)
+         call mpas_pool_get_array(thermalPool, 'iceOceanInterfaceTemperature', iceOceanInterfaceTemperature)
 
          if (nISMIP6OceanLayers > 0) then
             orig3dOceanMask(:,:) = 0
@@ -1459,7 +1465,8 @@ subroutine glc_import_mct(x2g_g, errorCode)
                call mpas_log_write("Unknown value for 'config_smb_source'", MPAS_LOG_ERR)
                errorCode = ior(errorCode, 1)
             endif
-            floatingBasalMassBal(i) = x2g_g % rAttr(index_x2g_Fogx_qiceli, n)
+            floatingBasalMassBal(i) = -(x2g_g % rAttr(index_x2g_Foxo_ismw, n) &
+                                      + x2g_g % rAttr(index_x2g_Foxo_frazil_li, n))
             if (nISMIP6OceanLayers > 0) then
                do iLev = 1, nISMIP6OceanLayers
                   fractionalMaskVal = x2g_g % rAttr(index_x2g_So_tf3d_mask(iLev), n)
@@ -1475,10 +1482,14 @@ subroutine glc_import_mct(x2g_g, errorCode)
                   endif
                enddo
             endif
-!            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)
-            !OceanDensity (i) = x2g_g % rAttr(, n)
+            basalOceanHeatFlux(i) = -x2g_g % rAttr(index_x2g_Fogx_qicehi, n)
+            fractionalMaskVal = x2g_g % rAttr(index_x2g_So_liflfrac, n)
+            if (fractionalMaskVal > config_min_floating_frac) then
+               iceOceanInterfaceTemperature(i) = x2g_g % rAttr(index_x2g_So_tfrz_isf, n) &
+                                                 / fractionalMaskVal
+            else
+               iceOceanInterfaceTemperature(i) = config_invalid_value_ocn_intr_temp
+            endif
          end do
 
          block => block % next
@@ -1516,13 +1527,12 @@ subroutine glc_export_mct(g2x_g, errorCode)
    integer, pointer :: nCellsSolve, nVertLevels
 
    real (kind=RKIND), dimension(:), pointer :: upperSurface
-   real (kind=RKIND), dimension(:), pointer :: layerThicknessFractions
    real (kind=RKIND), dimension(:), pointer :: thickness
+   real (kind=RKIND), dimension(:), pointer :: basalConductiveFlux
    real (kind=RKIND), dimension(:), pointer :: avgBareIceAblationApplied
    real (kind=RKIND), dimension(:), pointer :: avgCalvingFlux
    real (kind=RKIND), dimension(:), pointer :: avgFaceMeltFlux
    real (kind=RKIND), dimension(:), pointer :: avgFloatingBMBFlux
-   real (kind=RKIND), dimension(:,:), pointer :: temperature
    integer, dimension(:), pointer :: cellMask
    !-------------------------------------------------------------------
 
@@ -1545,9 +1555,8 @@ subroutine glc_export_mct(g2x_g, errorCode)
          call mpas_pool_get_subpool(block % structs, 'timeAveraging', timeAveragingPool)
 
          call mpas_pool_get_array(geometryPool, 'upperSurface', upperSurface)
-         call mpas_pool_get_array(meshPool, 'layerThicknessFractions', layerThicknessFractions)
          call mpas_pool_get_array(geometryPool, 'thickness', Thickness, timeLevel = 1)
-         call mpas_pool_get_array(thermalPool, 'temperature', temperature)
+         call mpas_pool_get_array(thermalPool, 'basalConductiveFlux', basalConductiveFlux)
 
          call mpas_pool_get_array(timeAveragingPool, 'avgBareIceAblationApplied', avgBareIceAblationApplied)
          call mpas_pool_get_array(timeAveragingPool, 'avgCalvingFlux', avgCalvingFlux)
@@ -1560,11 +1569,11 @@ subroutine glc_export_mct(g2x_g, errorCode)
             n = n + 1
 
             ! Fogg_rofl
-            g2x_g % rAttr(index_g2x_Fogg_rofl,n) = avgBareIceAblationApplied(i) 
+            g2x_g % rAttr(index_g2x_Fogg_rofl,n) = avgBareIceAblationApplied(i)
             ! Figg_rofi
             g2x_g % rAttr(index_g2x_Figg_rofi,n) = 0.0 ! placeholder
             ! Fogg_rofi
-            g2x_g % rAttr(index_g2x_Fogg_rofi,n) = avgCalvingFlux(i) 
+            g2x_g % rAttr(index_g2x_Fogg_rofi,n) = avgCalvingFlux(i)
             g2x_g % rAttr(index_g2x_Fogg_rofi,n) = g2x_g % rAttr(index_g2x_Fogg_rofi,n) + avgFaceMeltFlux(i)
             if (trim(config_basal_mass_bal_float) == 'ismip6') then
                ! if MALI is calculating ISMF, add that to rofi
@@ -1573,10 +1582,9 @@ subroutine glc_export_mct(g2x_g, errorCode)
                g2x_g % rAttr(index_g2x_Fogg_rofi,n) = g2x_g % rAttr(index_g2x_Fogg_rofi,n) - avgFloatingBMBFlux(i)
             endif
 
+            g2x_g % rAttr(index_g2x_Flgg_hflx, n) = basalConductiveFlux(i)
             g2x_g % rAttr(index_g2x_Sg_topo, n) = max(0.0, upperSurface(i)) !updated to avoid warning for values below sea level
-            g2x_g % rAttr(index_g2x_Sg_tbot, n) = temperature(nVertlevels,i) - SHR_CONST_TKTRIP
             ! layerThickness is not populated on init, so calculating like this instead:
-            g2x_g % rAttr(index_g2x_Sg_dztbot, n) = layerThicknessFractions(nVertLevels) * thickness(i) / 2.0
             g2x_g % rAttr(index_g2x_Sg_lithop, n) = Thickness(i)*SHR_CONST_RHOICE*SHR_CONST_G
 
             !!!Mask configurations
@@ -1648,27 +1656,27 @@ subroutine convert_seconds_to_timestamp(seconds, timeStamp)!{{{
       character (len=StrKIND), intent(out) :: timeStamp
       real (kind=RKIND) :: secondsPerHour, secondsPerMinute, remaining
       integer :: minutes, hours, secondsLeft
-   
+
       secondsPerHour = 3600
       secondsPerMinute = 60
-   
+
       if(seconds < 0 .or. seconds > 86400) then
         secondsLeft = 00
         minutes = 00
         hours = 00
       else
         hours = int(seconds/secondsPerHour)
         remaining = seconds - real(hours) * secondsPerHour
-   
+
         minutes = int(remaining/secondsPerMinute)
         remaining = remaining - real(minutes) * secondsPerMinute
-   
+
         secondsLeft = int(remaining)
       end if
-   
+
       write(timeStamp,"(a,i2.2,a,i2.2,a,i2.2)") "00_",hours,":",minutes,":",secondsLeft
       timeStamp = trim(timeStamp)
-   
+
    end subroutine convert_seconds_to_timestamp!}}}
 
    subroutine add_stream_attributes(stream_manager, domain)!{{{

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

@@ -1,5 +1,5 @@
 module glc_cpl_indices
-  
+
   use seq_flds_mod
   use mct_mod
 !  use glc_constants, only : glc_nec, glc_smb  ! TODO Will these be needed?  If so, need to add MPAS version
@@ -17,28 +17,20 @@ module glc_cpl_indices
   integer, public :: index_x2g_Flgl_qice = 0 !Ice sheet surface mass balance
   integer, public :: index_x2g_Sl_tsrf   = 0 !Ice sheet upper surface boundary temperature
 
-  integer, public :: index_x2g_So_blt    = 0 !Ice shelf boundary layer ocean temperature
-  integer, public :: index_x2g_So_bls    = 0 !Ice shelf boundary layer ocean salinity
-  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_tfrz_isf  = 0 !Ice shelf-ocean interface temperature
+  integer, public :: index_x2g_So_liflfrac  = 0 !Ice shelf fractional coverage from ocean
   integer, public :: index_x2g_So_tf3d(glc_nzoc_max)   = 0 !Ocean thermal forcing at specified z-levels
   integer, public :: index_x2g_So_tf3d_mask(glc_nzoc_max)   = 0 !mask of ocean thermal forcing at specified z-levels
-  integer, public :: index_x2g_Fogx_qiceli = 0 !Subshelf mass flux
   integer, public :: index_x2g_Fogx_qicehi = 0 !Subshelf heat flux for the ice sheet
+  integer, public :: index_x2g_Foxo_ismw = 0      !Ice shelf meltwater flux from ocean
+  integer, public :: index_x2g_Foxo_frazil_li = 0 !Frazil flux under ice shelves from ocean
 
   integer, public :: index_g2x_Fogg_rofi = 0 ! frozen runoff -> ocn
   integer, public :: index_g2x_Figg_rofi = 0 ! frozen runoff -> ice
   integer, public :: index_g2x_Fogg_rofl = 0 ! liquid runoff -> ocn
   integer, public :: index_g2x_Sg_topo = 0 ! surface elevation
   integer, public :: index_g2x_Flgg_hflx = 0 ! heat flux
-  integer, public :: index_g2x_Sg_tbot = 0 !bottom-layer temperature
-  integer, public :: index_g2x_Sg_dztbot = 0 !distance from bottom-layer temperature to ice-ocean interface (to compute T gradient)
   integer, public :: index_g2x_Sg_lithop = 0 !ice sheet lithostatic pressure
-  integer, public :: index_g2x_Fogx_qicelo = 0 !Subshelf liquid flux for ocean
-  integer, public :: index_g2x_Fogx_qiceho = 0 !Subshelf heat flux for the ocean
-  integer, public :: index_g2x_Sg_blis = 0 !Boundary layer interface salinity for ocean
-  integer, public :: index_g2x_Sg_blit = 0 !Boundary layer interface temperature for ocean
   integer, public :: index_g2x_Sg_icemask = 0 !complete grounded/floating ice mask
   integer, public :: index_g2x_Sg_ice_covered = 0 !ice covered ice mask
   integer, public :: index_g2x_Sg_icemask_grounded = 0 !grounded mask
@@ -75,9 +67,12 @@ subroutine glc_cpl_indices_set( )
     index_x2g_Sl_tsrf =                   mct_avect_indexra(x2g,'Sl_tsrf',perrwith='quiet')
     ! Surface temperature (deg C)
 
-    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')
+
+    ! Sub-ice-shelf fluxes from ocean
+    index_x2g_Foxo_ismw =                 mct_avect_indexra(x2g,'Foxo_ismw',perrwith='quiet')
+    index_x2g_Foxo_frazil_li =            mct_avect_indexra(x2g,'Foxo_frazil_li',perrwith='quiet')
+
 
     if (glc_nzoc > 0) then
        do iLev = 1, glc_nzoc
@@ -92,12 +87,6 @@ subroutine glc_cpl_indices_set( )
 
     !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)
-    index_x2g_So_blt =                    mct_avect_indexra(x2g,'So_blt',perrwith='quiet')
-    index_x2g_So_bls =                    mct_avect_indexra(x2g,'So_bls',perrwith='quiet')
-    index_x2g_So_htv =                    mct_avect_indexra(x2g,'So_htv',perrwith='quiet')
-    index_x2g_So_stv =                    mct_avect_indexra(x2g,'So_stv',perrwith='quiet')
-    index_g2x_Sg_tbot =                   mct_avect_indexra(g2x,'Sg_tbot',perrwith='quiet')
-    index_g2x_Sg_dztbot =                 mct_avect_indexra(g2x,'Sg_dztbot',perrwith='quiet')
     index_g2x_Sg_lithop =                 mct_avect_indexra(g2x,'Sg_lithop',perrwith='quiet')
 
     !Following block are GLC outputs for other components.
@@ -151,11 +140,9 @@ subroutine glc_cpl_indices_set( )
     ! but not evolving.  The GLC_TWO_WAY_COUPLING xml variable is what controls
     ! how this should be handled in CLM and the ice sheet model.
 
-    !Following block are subshelf melt routine outputs for ocean import
-    index_g2x_Fogx_qicelo = mct_avect_indexra(g2x,'Fogx_qicelo',perrwith='quiet')
-    index_g2x_Fogx_qiceho = mct_avect_indexra(g2x,'Fogx_qiceho',perrwith='quiet')
-    index_g2x_Sg_blis =     mct_avect_indexra(g2x,'Sg_blis',perrwith='quiet')
-    index_g2x_Sg_blit =     mct_avect_indexra(g2x,'Sg_blit',perrwith='quiet')
+    ! ice shelf-ocean interface temperature from the ocean
+    index_x2g_So_tfrz_isf = mct_avect_indexra(x2g,'So_tfrz_isf',perrwith='quiet')
+    index_x2g_So_liflfrac = mct_avect_indexra(x2g,'So_liflfrac',perrwith='quiet')
 
     call mct_aVect_clean(x2g)
     call mct_aVect_clean(g2x)

components/mpas-albany-landice/src/Registry.xml

@@ -490,11 +490,11 @@
 	   />
 		<nml_option name="config_TF_vertical_calc" type="character" default_value="seafloor" units="unitless"
 				  description="Method for calculating 2D thermal forcing from the 3d field. 'seafloor' uses the thermal forcing value immediately above the ocean floor (at grounding line depth); 'average' averages the full water column; '200-500' averages only between 200 - 500 m depth"
-		        possible_values="'seafloor', 'average', '200-500'"		  
+		        possible_values="'seafloor', 'average', '200-500'"
 		/>
 	   <nml_option name="config_TF_iceberg_melt" type="character" default_value="off" units="unitless"
 				  description="Option to modify ocean temperature in the upper portion of the water column (above depth specified by config_TF_iceberg_depth) to resemble iceberg melt-driven cooling. Done before calculating the 3d thermal forcing field if both config_ocean_data_extrapolation and config_calculate_thermal_forcing are also true. 'on' applies modification universally, while 'mask' uses binary mask 'icebergFjordMask' to apply the correction. "
-				  possible_values="'off', 'on', 'mask'"		  
+				  possible_values="'off', 'on', 'mask'"
 		/>
 		<nml_option name="config_TF_iceberg_depth" type="real" default_value="-175.0" units="m"
 				  description="Depth above which to iceberg meltwater is important. If config_TF_iceberg_melt is set to 'true' or 'mask', water above this depth will be altered to follow the submarine meltwater mixing line. Default value is taken from Hager et al. (2024)."
@@ -506,11 +506,19 @@
       <nml_option name="config_invalid_value_TF" type="real" default_value="1.0e36" units="unitless"
              description="value assigned to indicate invalid thermal forcing value when config_ocean_data_extrapolation is set to true"
              possible_values="Any real value"
-      />                     
+      />
 		<nml_option name="config_weight_value_cell" type="real" default_value="0.9" units="unitless"
 			    description="weight used to smooth horizontal extrapolation of ocean data field. Value close to 1 implies more weight of the current cell value versus the averaged value of the neighbouring cells around the cell"
 			    possible_values="any real value between 0 and 1"
 		/>
+		<nml_option name="config_invalid_value_ocn_intr_temp" type="real" default_value="1.0e36" units="unitless"
+				description="value assigned to indicate invalid values for the temperature at the ice-ocean interface"
+				possible_values="Any real value"
+		/>
+        <nml_option name="config_min_floating_frac" type="real" default_value="0.01" units="unitless"
+				description="The minimum floating fraction from the ocean for renormalizing the ice-ocean interface temperature"
+				possible_values="Any real value between 0 and 1"
+        />
 	</nml_record>
 
 	<nml_record name="physical_parameters" in_defaults="true">
@@ -1564,10 +1572,10 @@ is the value of that variable from the *previous* time level!
                 <var name="avgFloatingBMBFlux" type="real" dimensions="nCells Time" units="kg m^-2 s^-1"
                      description="time averaged floating basal mass balance"
                 />
-                <var name="avgBareIceAblation" type="real" dimensions="nCells Time" units="kg m^-2 s^-1" 
+                <var name="avgBareIceAblation" type="real" dimensions="nCells Time" units="kg m^-2 s^-1"
                      description="time averaged potential negative surface mass balance for melting of bare ice (using BareIceAblation from 'geometry' pool)"
                 />
-                <var name="avgBareIceAblationApplied" type="real" dimensions="nCells Time" units="kg m^-2 s^-1" 
+                <var name="avgBareIceAblationApplied" type="real" dimensions="nCells Time" units="kg m^-2 s^-1"
                      description="time averaged negative surface mass balance applied to melting of bare ice (using appliedBareIceAblation from 'geometry' pool)"
                 />
         </var_struct>
@@ -1772,12 +1780,18 @@ is the value of that variable from the *previous* time level!
                 <var name="basalHeatFlux" type="real" dimensions="nCells Time" units="W m^-2" default_value="0.06"
                      description="basal heat flux into the ice (positive upward)"
                 />
+                <var name="basalOceanHeatFlux" type="real" dimensions="nCells Time" units="W m^-2"
+                     description="ocean-driven basal heat flux under floating ice (positive downward)"
+                />
                 <var name="basalFrictionFlux" type="real" dimensions="nCells Time" units="W m^-2"
                      description="basal frictional heat flux into the ice (positive upward)"
                 />
                 <var name="heatDissipation" type="real" dimensions="nVertLevels nCells Time" units="C s^-1"
                      description="interior heat dissipation rate, divided by rhoi*cp_ice"
                 />
+                <var name="iceOceanInterfaceTemperature" type="real" dimensions="nCells Time" units="K"
+                     description="temperature at the ice shelf-ocean interface"
+                />
         </var_struct>
 
 <!-- ================ -->