Make depth at which to calc TF namelist-configurable

E3SM-Ocean-Discussion · May 13, 2024 · a91868e · a91868e
1 parent 46a76a3
commit a91868e
Show file tree

Hide file tree

Showing 7 changed files with 34 additions and 28 deletions.
diff --git a/components/mpas-albany-landice/driver/glc_comp_mct.F b/components/mpas-albany-landice/driver/glc_comp_mct.F
@@ -1410,7 +1410,7 @@ subroutine glc_import_mct(x2g_g, errorCode)
             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_tf300, 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)

diff --git a/components/mpas-albany-landice/driver/glc_cpl_indices.F b/components/mpas-albany-landice/driver/glc_cpl_indices.F
@@ -22,7 +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_tf300  = 0 !Ocean thermal forcing at 300m
+  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
 
@@ -71,7 +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_tf300  =                 mct_avect_indexra(x2g,'So_tf300',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)

diff --git a/components/mpas-ocean/driver/mpaso_cpl_indices.F b/components/mpas-ocean/driver/mpaso_cpl_indices.F
@@ -37,7 +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_tf300 !ocean thermal forcing at 300m
+  integer :: index_o2x_So_tf2d !ocean thermal forcing at predefined critical depth
 
 
   ! ocn -> drv  (BGC)
@@ -209,7 +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_tf300      = mct_avect_indexra(o2x,'So_tf300')
+    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')

diff --git a/components/mpas-ocean/driver/ocn_comp_mct.F b/components/mpas-ocean/driver/ocn_comp_mct.F
@@ -2703,7 +2703,7 @@ subroutine ocn_export_mct(o2x_o, errorCode)   !{{{
                                                avgRemovedIceRunoffFlux, &
                                                avgLandIceHeatFlux, &
                                                avgRemovedIceRunoffHeatFlux, &
-                                               avgThermalForcing300m
+                                               avgThermalForcingAtCritDepth
 
    real (kind=RKIND), dimension(:,:), pointer :: avgTracersSurfaceValue, avgSurfaceVelocity, &
                                                  avgSSHGradient, avgOceanSurfacePhytoC, &
@@ -2765,7 +2765,7 @@ subroutine ocn_export_mct(o2x_o, errorCode)   !{{{
      call mpas_pool_get_array(forcingPool, 'avgSurfaceVelocity', avgSurfaceVelocity)
      call mpas_pool_get_array(forcingPool, 'avgSSHGradient', avgSSHGradient)
      call mpas_pool_get_array(forcingPool, 'avgTotalFreshWaterTemperatureFlux', avgTotalFreshWaterTemperatureFlux)
-     call mpas_pool_get_array(forcingPool, 'avgThermalForcing300m', avgThermalForcing300m)
+     call mpas_pool_get_array(forcingPool, 'avgThermalForcingAtCritDepth', avgThermalForcingAtCritDepth)
 
      if ( frazilIceActive ) then
         call mpas_pool_get_array(forcingPool, 'seaIceEnergy', seaIceEnergy)
@@ -2942,7 +2942,7 @@ 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
-       o2x_o % rAttr(index_o2x_So_tf300, n) = avgThermalForcing300m(i)
+       o2x_o % rAttr(index_o2x_So_tf2d, n) = avgThermalForcingAtCritDepth(i)
 
 
        !Fyke: test

diff --git a/components/mpas-ocean/src/Registry.xml b/components/mpas-ocean/src/Registry.xml
@@ -770,6 +770,10 @@
 					description="If true, solid and liquid 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_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"
@@ -3896,8 +3900,8 @@
 			description="The time-averaged effective ocean density within ice shelves based on Archimedes' principle."
 			packages="landIceCouplingPKG"
 		/>
-		<var name="avgThermalForcing300m" type="real" dimensions="nCells Time" units="C"
-			description="The time-averaged thermal forcing at 300 m depth"
+		<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"

diff --git a/components/mpas-ocean/src/shared/mpas_ocn_time_average_coupled.F b/components/mpas-ocean/src/shared/mpas_ocn_time_average_coupled.F
@@ -60,7 +60,7 @@ subroutine ocn_time_average_coupled_init(forcingPool)!{{{
                                                     avgRemovedRiverRunoffFlux, avgRemovedIceRunoffFlux, &
                                                     avgLandIceHeatFlux, avgRemovedIceRunoffHeatFlux
 
-        real (kind=RKIND), dimension(:), pointer :: avgThermalForcing300m
+        real (kind=RKIND), dimension(:), pointer :: avgThermalForcingAtCritDepth
 
         integer :: iCell
         integer, pointer :: nAccumulatedCoupled, nCells
@@ -154,12 +154,12 @@ subroutine ocn_time_average_coupled_init(forcingPool)!{{{
            !$omp end parallel
         end if
 
-        call mpas_pool_get_array(forcingPool, 'avgThermalForcing300m', avgThermalForcing300m)
+        call mpas_pool_get_array(forcingPool, 'avgThermalForcingAtCritDepth', avgThermalForcingAtCritDepth)
 
         !$omp parallel
         !$omp do schedule(runtime)
         do iCell = 1, nCells
-           avgThermalForcing300m(iCell) = 0.0_RKIND
+           avgThermalForcingAtCritDepth(iCell) = 0.0_RKIND
         end do
         !$omp end do
         !$omp end parallel
@@ -272,13 +272,13 @@ subroutine ocn_time_average_coupled_accumulate(statePool, forcingPool, timeLevel
                                                     removedRiverRunoffFlux, avgRemovedRiverRunoffFlux, &
                                                     removedIceRunoffFlux, avgRemovedIceRunoffFlux, &
                                                     avgRemovedIceRunoffHeatFlux, &
-                                                    avgThermalForcing300m
+                                                    avgThermalForcingAtCritDepth
 
         type (mpas_pool_type), pointer :: tracersPool
 
         real (kind=RKIND), dimension(:,:,:), pointer :: activeTracers
         integer, pointer :: indexTemperature, indexSalinity
-        integer :: iLevel0300, iLevel
+        integer :: iLevelCritDepth, iLevel
         real (kind=RKIND) :: freezingTemp
 
         real (kind=RKIND), dimension(:,:,:), pointer :: &
@@ -417,17 +417,19 @@ subroutine ocn_time_average_coupled_accumulate(statePool, forcingPool, timeLevel
 
         end if
 
-        call mpas_pool_get_array(forcingPool, 'avgThermalForcing300m', avgThermalForcing300m)
+        call mpas_pool_get_array(forcingPool, 'avgThermalForcingAtCritDepth', avgThermalForcingAtCritDepth)
         call mpas_pool_get_subpool(statePool, 'tracers', tracersPool)
         call mpas_pool_get_array(tracersPool, 'activeTracers', activeTracers, 2)
         call mpas_pool_get_dimension(tracersPool, 'index_temperature', indexTemperature)
         call mpas_pool_get_dimension(tracersPool, 'index_salinity', indexSalinity)
 
-        ! find find vertical level that is just above the 300 m reference level
-        iLevel0300 = 1
+
+        ! find vertical level that is just above the critical depth reference level
+        ! this does not account for depression due to ice shelf cavities or sea ice
+        iLevelCritDepth = 1
         do iLevel = 1, nVertLevels
-           if(refBottomDepth(iLevel) > 300.0_RKIND) then
-              iLevel0300 = iLevel-1
+           if(refBottomDepth(iLevel) > config_2d_thermal_forcing_depth) then
+              iLevelCritDepth = iLevel-1
               exit
            end if
         end do
@@ -436,13 +438,13 @@ subroutine ocn_time_average_coupled_accumulate(statePool, forcingPool, timeLevel
         ! calculate thermal forcing at identified level for each cell
         do iCell = 1, nCells
            ! ignore cells that are too shallow
-           if (iLevel0300 <= maxLevelCell(iCell)) then
+           if (iLevelCritDepth <= maxLevelCell(iCell)) then
               ! this uses the level shallower than the reference level.  could interpolate instead
               ! note: assuming no LandIce cavity, but we may want to support that
-              freezingTemp = ocn_freezing_temperature(salinity=activeTracers(indexSalinity, iLevel0300, iCell), &
-                 pressure=pressure(iLevel0300, iCell),  inLandIceCavity=.false.)
-              avgThermalForcing300m(iCell) = ( avgThermalForcing300m(iCell) * nAccumulatedCoupled &
-                 + activeTracers(indexTemperature, iLevel0300, iCell) - freezingTemp ) / ( nAccumulatedCoupled + 1)
+              freezingTemp = ocn_freezing_temperature(salinity=activeTracers(indexSalinity, iLevelCritDepth, iCell), &
+                 pressure=pressure(iLevelCritDepth, iCell),  inLandIceCavity=.false.)
+              avgThermalForcingAtCritDepth(iCell) = ( avgThermalForcingAtCritDepth(iCell) * nAccumulatedCoupled &
+                 + activeTracers(indexTemperature, iLevelCritDepth, iCell) - freezingTemp ) / ( nAccumulatedCoupled + 1)
             end if
         end do
         !$omp end do

diff --git a/driver-mct/shr/seq_flds_mod.F90 b/driver-mct/shr/seq_flds_mod.F90
@@ -2987,12 +2987,12 @@ subroutine seq_flds_set(nmlfile, ID, infodata)
        attname  = 'So_rhoeff'
        call metadata_set(attname, longname, stdname, units)
 
-       name = 'So_tf300'
+       name = 'So_tf2d'
        call seq_flds_add(o2x_states,trim(name))
        call seq_flds_add(x2g_states,trim(name))
        call seq_flds_add(x2g_tf_states_from_ocn,trim(name))
-       longname = 'ocean thermal forcing at 300 m depth'
-       stdname  = 'ocean_thermal_forcing_at_300m'
+       longname = 'ocean thermal forcing at predefined critical depth'
+       stdname  = 'ocean_thermal_forcing_at_critical_depth'
        units    = 'C'
        attname  = name
        call metadata_set(attname, longname, stdname, units)