Adds brine rejection parameterization v1

code builds

components/mpas-ocean/src/shared/Makefile

@@ -149,7 +149,7 @@ mpas_ocn_tracer_hmix_redi.o: mpas_ocn_constants.o mpas_ocn_config.o
 
 mpas_ocn_high_freq_thickness_hmix_del2.o: mpas_ocn_constants.o mpas_ocn_config.o mpas_ocn_mesh.o
 
-mpas_ocn_tracer_nonlocalflux.o: mpas_ocn_constants.o mpas_ocn_config.o
+mpas_ocn_tracer_nonlocalflux.o: mpas_ocn_constants.o mpas_ocn_config.o mpas_ocn_diagnostics_variables.o mpas_ocn_mesh.o
 
 mpas_ocn_tracer_short_wave_absorption.o: mpas_ocn_tracer_short_wave_absorption_jerlov.o mpas_ocn_tracer_short_wave_absorption_variable.o mpas_ocn_constants.o mpas_ocn_config.o
 

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

@@ -668,22 +668,24 @@ subroutine ocn_surface_bulk_forcing_active_tracers(meshPool, forcingPool, tracer
                                                            - (snowFlux(iCell) + iceRunoffFlux(iCell)) &
                                                            * latent_heat_fusion_mks) * hflux_factor
 
-        ! Negative seaIceSalinityFlux is an extraction of salt from the ocean
-        ! So, we negate seaIceSalinityFlux when determining how much salt this flux needs.
-        requiredSalt = - seaIceSalinityFlux(iCell) * sflux_factor * dt / layerThickness(minLevelCell(iCell), iCell)
-        allowedSalt = min( 4.0_RKIND, tracerGroup(index_salinity_flux, minLevelCell(iCell), iCell) )
-
-        if ( allowedSalt < requiredSalt ) then
-           tracersSurfaceFluxRemoved(index_salinity_flux, iCell) = tracersSurfaceFluxRemoved(index_salinity_flux, iCell)  &
-                                                                 + ( 1 - allowedSalt / requiredSalt ) * seaIceSalinityFlux(iCell) &
-                                                                 * sflux_factor
-
-           tracersSurfaceFlux(index_salinity_flux, iCell) = tracersSurfaceFlux(index_salinity_flux, iCell)  &
-                                                          + ( allowedSalt / requiredSalt ) * seaIceSalinityFlux(iCell) &
-                                                          * sflux_factor
-        else
-           tracersSurfaceFlux(index_salinity_flux, iCell) = tracersSurfaceFlux(index_salinity_flux, iCell)  &
+        if (.not. config_use_brine_rejection_parameterization) then
+           ! Negative seaIceSalinityFlux is an extraction of salt from the ocean
+           ! So, we negate seaIceSalinityFlux when determining how much salt this flux needs.
+           requiredSalt = - seaIceSalinityFlux(iCell) * sflux_factor * dt / layerThickness(minLevelCell(iCell), iCell)
+           allowedSalt = min( 4.0_RKIND, tracerGroup(index_salinity_flux, minLevelCell(iCell), iCell) )
+
+           if ( allowedSalt < requiredSalt ) then
+              tracersSurfaceFluxRemoved(index_salinity_flux, iCell) = tracersSurfaceFluxRemoved(index_salinity_flux, iCell)  &
+                                                                    + ( 1 - allowedSalt / requiredSalt ) * seaIceSalinityFlux(iCell) &
+                                                                    * sflux_factor
+
+              tracersSurfaceFlux(index_salinity_flux, iCell) = tracersSurfaceFlux(index_salinity_flux, iCell)  &
+                                                             + ( allowedSalt / requiredSalt ) * seaIceSalinityFlux(iCell) &
+                                                             * sflux_factor
+           else
+              tracersSurfaceFlux(index_salinity_flux, iCell) = tracersSurfaceFlux(index_salinity_flux, iCell)  &
                                                           + seaIceSalinityFlux(iCell) * sflux_factor
+           end if
         end if
       end do
       !$omp end do

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

@@ -128,6 +128,7 @@ subroutine ocn_tend_thick(tendPool, forcingPool)!{{{
       ! pointers for the actual output variables within the pools above
       real (kind=RKIND), dimension(:), pointer, contiguous :: &
          surfaceThicknessFlux,     &! surface thickness flux
+         seaIceSalinityFlux,       &
          surfaceThicknessFluxRunoff, &! surface thickness flux from runoff
          surfaceThicknessFluxSubglacialRunoff ! surface thickness flux from runoff
 
@@ -627,6 +628,7 @@ subroutine ocn_tend_tracer(tendPool, statePool, forcingPool, &
 
       real (kind=RKIND), dimension(:), pointer, contiguous :: &
          penetrativeTemperatureFlux, &! heat flux penetrating below sfc
+         seaIceSalinityFlux, &
          tracerGroupExponentialDecayRate ! exp decay rate for forcing
 
       real (kind=RKIND), dimension(:,:), pointer, contiguous :: &
@@ -702,6 +704,8 @@ subroutine ocn_tend_tracer(tendPool, statePool, forcingPool, &
                                              penetrativeTemperatureFlux)
       call mpas_pool_get_array(forcingPool, 'fractionAbsorbed', &
                                              fractionAbsorbed)
+      call mpas_pool_get_array(forcingPool, 'seaIceSalinityFlux', &
+                                             seaIceSalinityFlux)
       call mpas_pool_get_array(forcingPool, 'fractionAbsorbedRunoff', &
                                              fractionAbsorbedRunoff)
       call mpas_pool_get_array(forcingPool, 'fractionAbsorbedSubglacialRunoff', &
@@ -1339,15 +1343,21 @@ subroutine ocn_tend_tracer(tendPool, statePool, forcingPool, &
                endif ! compute budgets
 
                if (isActiveTracer) then
-                  call ocn_tracer_nonlocalflux_tend(meshPool, &
+                  call ocn_tracer_nonlocalflux_tend( &
                                           vertNonLocalFlux, &
                                           nonLocalSurfaceTracerFlux, &
                                           tracerGroupTend, err)
                else ! not active
-                  call ocn_tracer_nonlocalflux_tend(meshPool, &
+                  call ocn_tracer_nonlocalflux_tend( &
                                           vertNonLocalFlux, &
                                           tracerGroupSurfaceFlux, &
                                           tracerGroupTend, err)
+                  if (config_use_brine_rejection_parameterization) then
+                     call ocn_brine_rejection_tendency(tracerGroupTend, &
+                                          layerThickness, &
+                                          seaIceSalinityFlux, &
+                                          indexSalinity, dt, err)
+                  end if
                endif ! active tracer
 
                if (computeBudgets) then

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

@@ -26,6 +26,8 @@ module ocn_tracer_nonlocalflux
    use mpas_pool_routines
    use ocn_constants
    use ocn_config
+   use ocn_diagnostics_variables
+   use ocn_mesh
 
    implicit none
    private
@@ -44,7 +46,8 @@ module ocn_tracer_nonlocalflux
    !--------------------------------------------------------------------
 
    public :: ocn_tracer_nonlocalflux_tend, &
-             ocn_tracer_nonlocalflux_init
+             ocn_tracer_nonlocalflux_init, &
+             ocn_brine_rejection_tendency
 
    !--------------------------------------------------------------------
    !
@@ -58,6 +61,84 @@ module ocn_tracer_nonlocalflux
 
 contains
 
+!***********************************************************************
+!
+!  routine ocn_brine_rejction_tend
+!
+!> \brief   Computes tendency term due to brine rejection
+!> \author  Luke Van Roekel
+!> \date    11/08/24
+!> \details
+!>  This routine computes the tendency for salinity from brine rejection, based on 
+!>      Nguyen et al 2009
+!
+!-----------------------------------------------------------------------
+
+   subroutine ocn_brine_rejection_tendency(tend, layerThickness, seaIceSalinityFlux, index_salinity, dt, err)
+      !----------------------------------------------------------------
+      !
+      ! variables
+      !
+      !----------------------------------------------------------------
+
+      real (kind=RKIND), dimension(:,:,:), intent(inout) :: &
+         tend !< tracer tendency, only will modify the salinity tend here
+
+      real (kind=RKIND), dimension(:), intent(in) :: &
+         seaIceSalinityFlux
+
+      real (kind=RKIND), dimension(:,:), intent(in) :: &
+         layerThickness
+
+      real, dimension(:), allocatable :: sis_flux
+
+      real (kind=RKIND) :: dt
+
+      integer, intent(in) :: index_salinity
+
+      integer :: err
+
+      !----------------------------------------------------------------
+      !
+      ! local variables
+      !
+      !----------------------------------------------------------------
+
+      integer :: iCell, k, nCells
+
+      real (kind=RKIND) :: remainder, z, integral, D, A, rejectedSalt
+
+      nCells = nCellsHalo(1)
+
+      allocate(sis_flux(nVertLevels+1))
+
+      do iCell = 1, nCells
+        if (seaIceSalinityFlux(iCell) > 0) then !only compute for salinity into ocean
+           sis_flux(:) = 0.0_RKIND
+           D = boundaryLayerDepth(iCell)
+           z = 0.0_RKIND
+           A = (config_brine_param_n + 1) / D**(config_brine_param_n+1)
+           k = minLevelCell(iCell)
+           ! need a salinity flux top and bottom and the flux into a cell is the divergence??
+           do while (z < D)
+              sis_flux(k) = A*z**config_brine_param_n*seaIceSalinityFlux(iCell)
+              z = z + layerThickness(k,iCell)
+              k = k + 1
+           end do
+           ! need to pick up the last bit that goes over
+           sis_flux(k) = 1.0_RKIND ! this should finish
+
+           do k=minLevelCell(iCell),maxLevelCell(iCell)
+              ! need to flip over so we take k+1 - k
+              tend(index_salinity, k, iCell) = tend(index_salinity, k, iCell) + &
+                 sis_flux(k+1) - sis_flux(k)
+           end do
+        end if
+      end do
+
+      deallocate(sis_flux)
+  end subroutine ocn_brine_rejection_tendency
+
 !***********************************************************************
 !
 !  routine ocn_tracer_nonlocalflux_tend
@@ -70,16 +151,13 @@ module ocn_tracer_nonlocalflux
 !
 !-----------------------------------------------------------------------
 
-   subroutine ocn_tracer_nonlocalflux_tend(meshPool, vertNonLocalFlux, surfaceTracerFlux, tend, err)!{{{
+   subroutine ocn_tracer_nonlocalflux_tend(vertNonLocalFlux, surfaceTracerFlux, tend, err)!{{{
       !-----------------------------------------------------------------
       !
       ! input variables
       !
       !-----------------------------------------------------------------
 
-      type (mpas_pool_type), intent(in) :: &
-         meshPool          !< Input: mesh information
-
       real (kind=RKIND), dimension(:,:), intent(in) :: &
         surfaceTracerFlux !< Input: surface tracer fluxes
 
@@ -110,9 +188,6 @@ subroutine ocn_tracer_nonlocalflux_tend(meshPool, vertNonLocalFlux, surfaceTrace
       !-----------------------------------------------------------------
 
       integer :: iCell, k, iTracer, nTracers, nCells
-      integer, pointer :: nVertLevels
-      integer, dimension(:), pointer :: nCellsArray
-      integer, dimension(:), pointer :: minLevelCell, maxLevelCell
       real (kind=RKIND) :: fluxTopOfCell, fluxBottomOfCell
 
       err = 0
@@ -121,14 +196,9 @@ subroutine ocn_tracer_nonlocalflux_tend(meshPool, vertNonLocalFlux, surfaceTrace
 
       call mpas_timer_start('non-local flux')
 
-      call mpas_pool_get_dimension(meshPool, 'nCellsArray', nCellsArray)
-      call mpas_pool_get_dimension(meshPool, 'nVertLevels', nVertLevels)
       nTracers = size(tend, dim=1)
 
-      call mpas_pool_get_array(meshPool, 'minLevelCell', minLevelCell)
-      call mpas_pool_get_array(meshPool, 'maxLevelCell', maxLevelCell)
-
-      nCells = nCellsArray( 1 )
+      nCells = nCellsHalo( 1 )
 
       !$omp parallel
       !$omp do schedule(runtime) private(k, iTracer, fluxTopOfCell, fluxBottomOfCell)