Adds new ocean to ice coupling fields for frazil ice formation

components/mpas-ocean/driver/ocn_comp_mct.F

@@ -3166,6 +3166,7 @@ subroutine ocn_export_mct(o2x_o, errorCode)   !{{{
    integer, dimension(:), pointer :: landIceMask
 
    real (kind=RKIND), dimension(:), pointer :: seaIceEnergy, accumulatedFrazilIceMass, frazilSurfacePressure, &
+                                               accumulatedFrazilIceEnergy, &
                                                avgTotalFreshWaterTemperatureFlux, &
                                                avgCO2_gas_flux, DMSFlux, surfaceUpwardCO2Flux, &
                                                avgOceanSurfaceDIC, &
@@ -3286,6 +3287,7 @@ subroutine ocn_export_mct(o2x_o, errorCode)   !{{{
         call mpas_pool_get_array(forcingPool, 'seaIceEnergy', seaIceEnergy)
         call mpas_pool_get_array(forcingPool, 'frazilSurfacePressure', frazilSurfacePressure)
         call mpas_pool_get_array(statePool, 'accumulatedFrazilIceMass', accumulatedFrazilIceMass, 1)
+        call mpas_pool_get_array(statePool, 'accumulatedFrazilIceEnergy', accumulatedFrazilIceEnergy, 1)
      end if
 
      ! Cryo fields
@@ -3405,8 +3407,10 @@ subroutine ocn_export_mct(o2x_o, errorCode)   !{{{
 
              o2x_o(index_o2x_Fioo_q, n)  = seaIceEnergy(i) / ocn_cpl_dt
              o2x_o(index_o2x_Fioo_frazil, n) = accumulatedFrazilIceMass(i) / ocn_cpl_dt
-             o2x_o(index_o2x_Fioo_frazils, n) = 4.0_RKIND
-             o2x_o(index_o2x_Fioo_frazilh, n) = 0.0_RKIND
+             o2x_o(index_o2x_Fioo_frazils, n) = accumulatedFrazilIceMass(i)/ 1000.0_RKIND * &
+                config_frazil_sea_ice_reference_salinity / ocn_cpl_dt
+             !defined positive to be consistent with Fioo_q when frazil is formed
+             o2x_o(index_o2x_Fioo_frazilh, n) = accumulatedFrazilIceEnergy(i)/ ocn_cpl_dt 
 
           else
 
@@ -3424,6 +3428,7 @@ subroutine ocn_export_mct(o2x_o, errorCode)   !{{{
           ! Reset SeaIce Energy and Accumulated Frazil Ice
           seaIceEnergy(i) = 0.0_RKIND
           accumulatedFrazilIceMass(i) = 0.0_RKIND
+          accumulatedFrazilIceEnergy(i) = 0.0_RKIND
           frazilSurfacePressure(i) = 0.0_RKIND
        end if
 

components/mpas-ocean/src/Registry.xml

@@ -2146,6 +2146,7 @@
 			<var_struct name="ecosysAuxiliary"/>
 			<var name="accumulatedFrazilIceMass"/>
 			<var name="accumulatedFrazilIceSalinity"/>
+			<var name="accumulatedFrazilIceEnergy"/>
 			<var name="accumulatedLandIceMass"/>
 			<var name="accumulatedLandIceHeat"/>
 			<var name="accumulatedLandIceFrazilMass"/>
@@ -2517,6 +2518,10 @@
 			 description="Mass per unit area of frazil ice produced. Reset to zero at each coupling interval"
 			 packages="frazilIce"
 		/>
+		<var name="accumulatedFrazilIceEnergy" type="real" dimensions="nCells Time" units="J m^-2"
+			 description="Accumulated frazil ice energy per unit area. Reset to zero at each coupling interval"
+			 packages="frazilIce"
+		/>
 		<var name="accumulatedFrazilIceSalinity" type="real" dimensions="nCells Time" units="kg m^-2"
 			 description="Salinity associated with accumulatedFrazilIceMass. Reset to zero at each coupling interval"
 			 packages="frazilIce"

components/mpas-ocean/src/analysis_members/Registry_conservation_check.xml

@@ -86,7 +86,7 @@
 			description="Iceberg heat flux from coupler, integrated in space and time. Positive into the ocean."
 		/>
 		<var name="accumulatedFrazilHeatFlux" type="real" dimensions="Time" units="W m^-2"
-			description="Heat flux from frazil, integrated in space and time, given to couipler. Positive into the ocean."
+			description="Heat flux from frazil, integrated in space and time, given to coupler. Positive into the ocean."
 		/>
 		<var name="accumulatedLandIceHeatFlux" type="real" dimensions="Time" units="W m^-2"
 			description="Land ice heat flux, integrated in space and time. Positive into the ocean."

components/mpas-ocean/src/analysis_members/mpas_ocn_conservation_check.F

@@ -462,6 +462,8 @@ subroutine energy_conservation(domain, err)
            icebergTemperatureFlux, &
            accumulatedFrazilIceMassNew, &
            accumulatedFrazilIceMassOld, &
+           accumulatedFrazilIceEnergyNew, &
+           accumulatedFrazilIceEnergyOld, &
            accumulatedLandIceFrazilMassNew, &
            accumulatedLandIceFrazilMassOld
 
@@ -565,6 +567,8 @@ subroutine energy_conservation(domain, err)
             call mpas_pool_get_array(tracersSurfaceFluxPool, 'activeTracersSurfaceFluxSubglacialRunoff', activeTracersSurfaceFluxSubglacialRunoff)
             call mpas_pool_get_array(statePool, 'accumulatedFrazilIceMass', accumulatedFrazilIceMassNew, 2)
             call mpas_pool_get_array(statePool, 'accumulatedFrazilIceMass', accumulatedFrazilIceMassOld, 1)
+            call mpas_pool_get_array(statePool, 'accumulatedFrazilIceEnergy', accumulatedFrazilIceEnergyNew, 2)
+            call mpas_pool_get_array(statePool, 'accumulatedFrazilIceEnergy', accumulatedFrazilIceEnergyOld, 1)
 
             do iCell = 1, nCellsSolve
 
@@ -599,8 +603,10 @@ subroutine energy_conservation(domain, err)
                do iCell = 1, nCellsSolve
                   ! Frazil ice mass is negative. Negative coefficient makes heat
                   ! flux positive, because freezing ice releases heat.
-                  sumArray(10) = sumArray(10) - areaCell(iCell) * config_frazil_heat_of_fusion &
-                                 * (accumulatedFrazilIceMassNew(iCell) - accumulatedFrazilIceMassOld(iCell))/dt
+                  !sumArray(10) = sumArray(10) - areaCell(iCell) * config_frazil_heat_of_fusion &
+                  !               * (accumulatedFrazilIceMassNew(iCell) - accumulatedFrazilIceMassOld(iCell))/dt
+                  sumArray(10) = sumArray(10) - areaCell(iCell) &
+                                  * (accumulatedFrazilIceEnergyNew(iCell) - accumulatedFrazilIceEnergyOld(iCell))/dt                        
                enddo
             end if
 
@@ -749,7 +755,8 @@ subroutine energy_conservation(domain, err)
             s = 0.0_RKIND
             call mpas_log_write('HEAT FLUXES: explicit')
             call mpas_log_write('MPAS-Ocean name            W (raw sum)    coupler name    short name         W/m^2 (flux/A)')
-v=accumulatedFrazilHeatFlux          ; write(m,"('frazilMass * LHF /dt     ',es16.8,' o2x_Fioo_q      hfreeze         ',f16.8)") v,v/A; call mpas_log_write(m); s=s+v
+!v=accumulatedFrazilHeatFlux          ; write(m,"('frazilMass * LHF /dt     ',es16.8,' o2x_Fioo_q      hfreeze         ',f16.8)") v,v/A; call mpas_log_write(m); s=s+v
+v=accumulatedFrazilHeatFlux          ; write(m,"('frazilMass * LHF /dt     ',es16.8,' o2x_Fioo_frazilh hfreeze         ',f16.8)") v,v/A; call mpas_log_write(m); s=s+v
 v=accumulatedSeaIceHeatFlux          ; write(m,"('seaIceHeatFlux           ',es16.8,' x2o_Fioi_melth  hmelt           ',f16.8)") v,v/A; call mpas_log_write(m); s=s+v
 v=accumulatedShortWaveHeatFlux       ; write(m,"('shortWaveHeatFlux        ',es16.8,' x2o_Foxx_swnet  hnetsw          ',f16.8)") v,v/A; call mpas_log_write(m); s=s+v
 v=accumulatedLongWaveHeatFluxDown    ; write(m,"('longWaveHeatFluxDown     ',es16.8,' x2o_Faxa_lwdn   hlwdn           ',f16.8)") v,v/A; call mpas_log_write(m); s=s+v

components/mpas-ocean/src/shared/mpas_ocn_frazil_forcing.F

@@ -462,10 +462,13 @@ subroutine ocn_frazil_forcing_build_arrays(domain, meshPool, forcingPool, stateP
       real (kind=RKIND) :: newFrazilIceThickness, newThicknessWeightedSaltContent
       real (kind=RKIND) :: sumNewFrazilIceThickness, sumNewThicknessWeightedSaltContent
       real (kind=RKIND) :: meltedFrazilIceThickness, meltedThicknessWeightedSaltContent
+      real (kind=RKIND) :: sumNewFrazilEnergy, newFrazilEnergy, meltedFrazilEnergy
       real (kind=RKIND) :: oceanFreezingTemperature, frazilSalinity
 
       real (kind=RKIND), pointer, dimension(:)     :: accumulatedFrazilIceMassNew
       real (kind=RKIND), pointer, dimension(:)     :: accumulatedFrazilIceMassOld
+      real (kind=RKIND), pointer, dimension(:)     :: accumulatedFrazilIceEnergyNew
+      real (kind=RKIND), pointer, dimension(:)     :: accumulatedFrazilIceEnergyOld
       real (kind=RKIND), pointer, dimension(:)     :: accumulatedFrazilIceSalinityNew
       real (kind=RKIND), pointer, dimension(:)     :: accumulatedFrazilIceSalinityOld
       real (kind=RKIND), pointer, dimension(:)     :: accumulatedLandIceFrazilMassNew
@@ -515,6 +518,8 @@ subroutine ocn_frazil_forcing_build_arrays(domain, meshPool, forcingPool, stateP
       call mpas_pool_get_array(statePool, 'layerThickness', layerThickness, 1)
       call mpas_pool_get_array(statePool, 'accumulatedFrazilIceMass', accumulatedFrazilIceMassNew, 2)
       call mpas_pool_get_array(statePool, 'accumulatedFrazilIceMass', accumulatedFrazilIceMassOld, 1)
+      call mpas_pool_get_array(statePool, 'accumulatedFrazilIceEnergy', accumulatedFrazilIceEnergyNew, 2)
+      call mpas_pool_get_array(statePool, 'accumulatedFrazilIceEnergy', accumulatedFrazilIceEnergyOld, 1)
       call mpas_pool_get_array(statePool, 'accumulatedFrazilIceSalinity', accumulatedFrazilIceSalinityNew, 2)
       call mpas_pool_get_array(statePool, 'accumulatedFrazilIceSalinity', accumulatedFrazilIceSalinityOld, 1)
       call mpas_pool_get_array(statePool, 'accumulatedLandIceFrazilMass', accumulatedLandIceFrazilMassNew, 2)
@@ -552,7 +557,8 @@ subroutine ocn_frazil_forcing_build_arrays(domain, meshPool, forcingPool, stateP
       !$omp         sumNewThicknessWeightedSaltContent, oceanFreezingTemperature, potential, &
       !$omp         freezingEnergy, meltingEnergy, frazilSalinity, newFrazilIceThickness, &
       !$omp         newThicknessWeightedSaltContent, meltedFrazilIceThickness, &
-      !$omp         meltedThicknessWeightedSaltContent)
+      !$omp         meltedThicknessWeightedSaltContent, newFrazilEnergy, meltedFrazilEnergy, &
+      !$omp         sumNewFrazilEnergy, newFrazilEnergy)
       do iCell=1,nCells
 
         underLandIce = .false.
@@ -604,6 +610,8 @@ subroutine ocn_frazil_forcing_build_arrays(domain, meshPool, forcingPool, stateP
         ! initialize the sum of new frazil ice created
         sumNewFrazilIceThickness = 0.0_RKIND
         sumNewThicknessWeightedSaltContent = 0.0_RKIND
+        sumNewFrazilEnergy = 0.0_RKIND
+        newFrazilEnergy = 0.0_RKIND
 
         ! loop from maximum depth of frazil creation to surface
         do k = kBottomFrazil, minLevelCell(iCell), -1
@@ -632,31 +640,29 @@ subroutine ocn_frazil_forcing_build_arrays(domain, meshPool, forcingPool, stateP
 
             ! new frazil ice formation measured in meters
             newFrazilIceThickness = freezingEnergy / (config_frazil_heat_of_fusion * config_frazil_ice_density)
-
             ! limit the frazil formed appropriately
             newFrazilIceThickness = min(newFrazilIceThickness, layerThickness(k,iCell) * config_frazil_fractional_thickness_limit)
 
             ! Determine new salt in frazil
             newThicknessWeightedSaltContent = newFrazilIceThickness * frazilSalinity
-
             ! compute tendency to thickness, temperature and salinity
             ! layerTendency is scaled so that mass of ice created == mass of ocean water removed
 
             ! layer thickness decreased due to creation of frazil
             ! note: -- this has to be density (not rho_sw) to keep buoyancy equal
             frazilLayerThicknessTendency(k,iCell) = - newFrazilIceThickness * config_frazil_ice_density / density(k,iCell) / dt
-
             ! salt is extracted with the frazil
             frazilSalinityTendency(k,iCell) = - newThicknessWeightedSaltContent / dt
 
             ! ocean fluid temperature is warmed due to creation of frazil
             frazilTemperatureTendency(k,iCell) =  + ( newFrazilIceThickness * config_frazil_heat_of_fusion &
                                                * config_frazil_ice_density ) / (config_specific_heat_sea_water * rho_sw) / dt
-
-            ! keep track of sum of frazil ice
+            ! Store the latent energy of the new frazil 
+            newFrazilEnergy = newFrazilIceThickness * config_frazil_ice_density * config_frazil_heat_of_fusion
+            !keep track of sum of frazil ice
             sumNewFrazilIceThickness = sumNewFrazilIceThickness + newFrazilIceThickness
             sumNewThicknessWeightedSaltContent = sumNewThicknessWeightedSaltContent + newThicknessWeightedSaltContent
-
+            sumNewFrazilEnergy = sumNewFrazilEnergy + newFrazilEnergy
           else
             ! ocean water is warm enough to melt frazil
 
@@ -688,18 +694,24 @@ subroutine ocn_frazil_forcing_build_arrays(domain, meshPool, forcingPool, stateP
               frazilSalinityTendency(k,iCell) = + meltedThicknessWeightedSaltContent/ dt
 
               ! ocean fluid temperature is cooled due to melting of frazil
-              frazilTemperatureTendency(k,iCell) = - ( meltedFrazilIceThickness * config_frazil_heat_of_fusion &
+              frazilTemperatureTendency(k,iCell) =  - ( meltedFrazilIceThickness * config_frazil_heat_of_fusion &
                                                  * config_frazil_ice_density ) / (config_specific_heat_sea_water * rho_sw) &
                                                  / dt
 
               ! ocean fluid temperature is modified due to mixing of existing water with meltwater at the freezing point
+              ! more consistent with the layer tendency
+              !frazilTemperatureTendency(k,iCell) = frazilTemperatureTendency(k,iCell) + &
+              !     meltedFrazilIceThickness * oceanFreezingTemperature * config_frazil_ice_density / density(k,iCell) / dt
+              ! original formulation
               frazilTemperatureTendency(k,iCell) = frazilTemperatureTendency(k,iCell) + &
                    meltedFrazilIceThickness * oceanFreezingTemperature / dt
-
               ! keep track of new frazil ice
               sumNewThicknessWeightedSaltContent = max(0.0_RKIND,sumNewThicknessWeightedSaltContent - meltedThicknessWeightedSaltContent)
               sumNewFrazilIceThickness = max(0.0_RKIND,sumNewFrazilIceThickness - meltedFrazilIceThickness)
 
+              meltedFrazilEnergy = meltedFrazilIceThickness * config_frazil_ice_density * config_frazil_heat_of_fusion
+
+              sumNewFrazilEnergy = max(0.0_RKIND,sumNewFrazilEnergy - meltedFrazilEnergy)
             endif  ! if (sumNewFrazilIceThickness > 0.0_RKIND)
 
           endif  ! if (freezingEnergy < 0)
@@ -711,6 +723,9 @@ subroutine ocn_frazil_forcing_build_arrays(domain, meshPool, forcingPool, stateP
         accumulatedFrazilIceMassNew(iCell) = accumulatedFrazilIceMassOld(iCell) + sumNewFrazilIceThickness &
            * config_frazil_ice_density
 
+        !frazilIceEnergy > 0
+        accumulatedFrazilIceEnergyNew(iCell) = accumulatedFrazilIceEnergyOld(iCell) + sumNewFrazilEnergy
+
         accumulatedFrazilIceSalinityNew(iCell) = accumulatedFrazilIceSalinityOld(iCell) + &
            (1-usePorosity) * sumNewThicknessWeightedSaltContent
 

components/mpas-seaice/bld/namelist_files/namelist_definition_mpassi.xml

@@ -533,7 +533,7 @@ Default: Defined in namelist_defaults.xml
 	category="forcing" group="forcing">
 Type of frazil ice coupling.
 
-Valid values: 'external', 'internal', 'fresh_ice_correction'
+Valid values: 'external', 'internal', 'fresh_ice_correction', 'omega-fluxes'
 Default: Defined in namelist_defaults.xml
 </entry>
 

components/mpas-seaice/cime_config/buildnml

@@ -852,6 +852,9 @@ def buildnml(case, caseroot, compname):
             lines.append('    <var name="snowVolumeCell"/>')
             lines.append('    <var name="uVelocityGeo"/>')
             lines.append('    <var name="vVelocityGeo"/>')
+            lines.append('    <var name="frazilMassFlux"/>')
+            lines.append('    <var name="frazilSaltFlux"/>')
+            lines.append('    <var name="frazilEnthalpyFlux"/>')
             if ("WW3" in compset or "DWAV%FSD" in compset or "DWAV%ZEROS" in compset):
                 lines.append('    <var name="floeSizeDiameter"/> ')
                 lines.append('    <var name="waveSpectra"/> ')
@@ -913,6 +916,9 @@ def buildnml(case, caseroot, compname):
             lines.append('    <var name="airTemperature"/>')
             lines.append('    <var name="congelation"/>')
             lines.append('    <var name="frazilFormation"/>')
+            lines.append('    <var name="frazilMassFlux"/>')
+            lines.append('    <var name="frazilSaltFlux"/>')
+            lines.append('    <var name="frazilEnthalpyFlux"/>')
             lines.append('    <var name="snowiceFormation"/>')
             lines.append('    <var name="snowMelt"/>')
             lines.append('    <var name="surfaceIceMelt"/>')