Skip to content

Commit

Permalink
Brine plume (#401)
Browse files Browse the repository at this point in the history 
* Salt data structures

* First steps at brine plume: pass info from SIS2

* The brine plume parameterization,

- including now passing the dimensional scaling tests.

* Fix problem when running Tidal_bay case with gnu.

* Avoiding visc_rem issues inside land mask.

Tweaking the brine plume code.

* Using the proper MLD in the brine plumes

- it now works better on restart

* Always including MLD in call to applyBoundary...

- I could move it up and make it not optional.

* Adding some OpenMP directives to brine plumes
  • Loading branch information
@kshedstrom
kshedstrom authored Jul 25, 2023
1 parent 2f6e86e commit 147ddf1
Show file tree
Hide file tree
Showing 7 changed files with 132 additions and 38 deletions.
13 changes: 12 additions & 1 deletion config_src/drivers/FMS_cap/MOM_surface_forcing_gfdl.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -62,7 +62,7 @@ module MOM_surface_forcing_gfdl
!! from MOM_domains) to indicate the staggering of
!! the winds that are being provided in calls to
!! update_ocean_model.
logical :: use_temperature !< If true, temp and saln used as state variables
logical :: use_temperature !< If true, temp and saln used as state variables.
real :: wind_stress_multiplier !< A multiplier applied to incoming wind stress [nondim].

real :: Rho0 !< Boussinesq reference density [R ~> kg m-3]
Expand Down Expand Up @@ -175,6 +175,7 @@ module MOM_surface_forcing_gfdl
real, pointer, dimension(:,:) :: t_flux =>NULL() !< sensible heat flux [W m-2]
real, pointer, dimension(:,:) :: q_flux =>NULL() !< specific humidity flux [kg m-2 s-1]
real, pointer, dimension(:,:) :: salt_flux =>NULL() !< salt flux [kg m-2 s-1]
real, pointer, dimension(:,:) :: excess_salt =>NULL() !< salt left behind by brine rejection [kg m-2 s-1]
real, pointer, dimension(:,:) :: lw_flux =>NULL() !< long wave radiation [W m-2]
real, pointer, dimension(:,:) :: sw_flux_vis_dir =>NULL() !< direct visible sw radiation [W m-2]
real, pointer, dimension(:,:) :: sw_flux_vis_dif =>NULL() !< diffuse visible sw radiation [W m-2]
Expand Down Expand Up @@ -304,6 +305,8 @@ subroutine convert_IOB_to_fluxes(IOB, fluxes, index_bounds, Time, valid_time, G,
call safe_alloc_ptr(fluxes%heat_added,isd,ied,jsd,jed)
call safe_alloc_ptr(fluxes%salt_flux_added,isd,ied,jsd,jed)

if (associated(IOB%excess_salt)) call safe_alloc_ptr(fluxes%salt_left_behind,isd,ied,jsd,jed)

do j=js-2,je+2 ; do i=is-2,ie+2
fluxes%TKE_tidal(i,j) = CS%TKE_tidal(i,j)
fluxes%ustar_tidal(i,j) = CS%ustar_tidal(i,j)
Expand Down Expand Up @@ -576,6 +579,11 @@ subroutine convert_IOB_to_fluxes(IOB, fluxes, index_bounds, Time, valid_time, G,
call check_mask_val_consistency(IOB%salt_flux(i-i0,j-j0), G%mask2dT(i,j), i, j, 'salt_flux', G)
enddo ; enddo
endif
if (associated(IOB%excess_salt)) then
do j=js,je ; do i=is,ie
fluxes%salt_left_behind(i,j) = G%mask2dT(i,j)*(kg_m2_s_conversion*IOB%excess_salt(i-i0,j-j0))
enddo ; enddo
endif

!#CTRL# if (associated(CS%ctrl_forcing_CSp)) then
!#CTRL# do j=js,je ; do i=is,ie
Expand Down Expand Up @@ -1729,6 +1737,9 @@ subroutine ice_ocn_bnd_type_chksum(id, timestep, iobt)
if (associated(iobt%mass_berg)) then
chks = field_chksum( iobt%mass_berg ) ; if (root) write(outunit,100) 'iobt%mass_berg ', chks
endif
if (associated(iobt%excess_salt)) then
chks = field_chksum( iobt%excess_salt ) ; if (root) write(outunit,100) 'iobt%excess_salt ', chks
endif
100 FORMAT(" CHECKSUM::",A20," = ",Z20)

call coupler_type_write_chksums(iobt%fluxes, outunit, 'iobt%')
Expand Down
9 changes: 9 additions & 0 deletions docs/zotero.bib
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2738,3 +2738,12 @@ @article{kraus1967
journal = {Tellus}
}

@article{Nguyen2009,
doi = {10.1029/2008JC005121},
year = {2009},
journal = {JGR Oceans},
volume = {114},
author = {A. T. Nguyen and D. Menemenlis and R. Kwok},
title = {Improved modeling of the Arctic halocline with a subgrid-scale brine rejection parameterization},
pages = {C11014}
}
8 changes: 4 additions & 4 deletions src/core/MOM_continuity_PPM.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -378,9 +378,9 @@ subroutine zonal_mass_flux(u, h_in, uh, dt, G, GV, US, CS, LB, OBC, por_face_are
dx_E = ratio_max(G%areaT(i+1,j), G%dy_Cu(I,j), 1000.0*G%dxT(i+1,j))
else ; dx_W = G%dxT(i,j) ; dx_E = G%dxT(i+1,j) ; endif

if (du_max_CFL(I) * visc_rem(I,k) > dx_W*CFL_dt - u(I,j,k)) &
if (du_max_CFL(I) * visc_rem(I,k) > dx_W*CFL_dt - u(I,j,k)*G%mask2dCu(I,j)) &
du_max_CFL(I) = (dx_W*CFL_dt - u(I,j,k)) / visc_rem(I,k)
if (du_min_CFL(I) * visc_rem(I,k) < -dx_E*CFL_dt - u(I,j,k)) &
if (du_min_CFL(I) * visc_rem(I,k) < -dx_E*CFL_dt - u(I,j,k)*G%mask2dCu(I,j)) &
du_min_CFL(I) = -(dx_E*CFL_dt + u(I,j,k)) / visc_rem(I,k)
enddo ; enddo
endif
Expand Down Expand Up @@ -1201,9 +1201,9 @@ subroutine meridional_mass_flux(v, h_in, vh, dt, G, GV, US, CS, LB, OBC, por_fac
dy_S = ratio_max(G%areaT(i,j), G%dx_Cv(I,j), 1000.0*G%dyT(i,j))
dy_N = ratio_max(G%areaT(i,j+1), G%dx_Cv(I,j), 1000.0*G%dyT(i,j+1))
else ; dy_S = G%dyT(i,j) ; dy_N = G%dyT(i,j+1) ; endif
if (dv_max_CFL(i) * visc_rem(i,k) > dy_S*CFL_dt - v(i,J,k)) &
if (dv_max_CFL(i) * visc_rem(i,k) > dy_S*CFL_dt - v(i,J,k)*G%mask2dCv(i,J)) &
dv_max_CFL(i) = (dy_S*CFL_dt - v(i,J,k)) / visc_rem(i,k)
if (dv_min_CFL(i) * visc_rem(i,k) < -dy_N*CFL_dt - v(i,J,k)) &
if (dv_min_CFL(i) * visc_rem(i,k) < -dy_N*CFL_dt - v(i,J,k)*G%mask2dCv(i,J)) &
dv_min_CFL(i) = -(dy_N*CFL_dt + v(i,J,k)) / visc_rem(i,k)
enddo ; enddo
endif
Expand Down
16 changes: 11 additions & 5 deletions src/core/MOM_forcing_type.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -135,8 +135,10 @@ module MOM_forcing_type
real, pointer, dimension(:,:) :: &
salt_flux => NULL(), & !< net salt flux into the ocean [R Z T-1 ~> kgSalt m-2 s-1]
salt_flux_in => NULL(), & !< salt flux provided to the ocean from coupler [R Z T-1 ~> kgSalt m-2 s-1]
salt_flux_added => NULL() !< additional salt flux from restoring or flux adjustment before adjustment
salt_flux_added => NULL(), & !< additional salt flux from restoring or flux adjustment before adjustment
!! to net zero [R Z T-1 ~> kgSalt m-2 s-1]
salt_left_behind => NULL() !< salt left in ocean at the surface from brine rejection
!! [R Z T-1 ~> kgSalt m-2 s-1]

! applied surface pressure from other component models (e.g., atmos, sea ice, land ice)
real, pointer, dimension(:,:) :: p_surf_full => NULL()
Expand Down Expand Up @@ -746,15 +748,15 @@ subroutine extractFluxes1d(G, GV, US, fluxes, optics, nsw, j, dt, &
endif

! Salt fluxes
Net_salt(i) = 0.0
if (do_NSR) Net_salt_rate(i) = 0.0
net_salt(i) = 0.0
if (do_NSR) net_salt_rate(i) = 0.0
! Convert salt_flux from kg (salt)/(m^2 * s) to
! Boussinesq: (ppt * m)
! non-Bouss: (g/m^2)
if (associated(fluxes%salt_flux)) then
Net_salt(i) = (scale * dt * (1000.0*US%ppt_to_S * fluxes%salt_flux(i,j))) * GV%RZ_to_H
net_salt(i) = (scale * dt * (1000.0*US%ppt_to_S * fluxes%salt_flux(i,j))) * GV%RZ_to_H
!Repeat above code for 'rate' term
if (do_NSR) Net_salt_rate(i) = (scale * 1. * (1000.0*US%ppt_to_S * fluxes%salt_flux(i,j))) * GV%RZ_to_H
if (do_NSR) net_salt_rate(i) = (scale * 1. * (1000.0*US%ppt_to_S * fluxes%salt_flux(i,j))) * GV%RZ_to_H
endif

! Diagnostics follow...
Expand Down Expand Up @@ -1947,6 +1949,10 @@ subroutine register_forcing_type_diags(Time, diag, US, use_temperature, handles,
diag%axesT1,Time,'Salt flux into ocean at surface due to restoring or flux adjustment', &
units='kg m-2 s-1', conversion=US%RZ_T_to_kg_m2s)

handles%id_saltFluxAdded = register_diag_field('ocean_model', 'salt_left_behind', &
diag%axesT1,Time,'Salt left in ocean at surface due to ice formation', &
units='kg m-2 s-1', conversion=US%RZ_T_to_kg_m2s)

handles%id_saltFluxGlobalAdj = register_scalar_field('ocean_model', &
'salt_flux_global_restoring_adjustment', Time, diag, &
'Adjustment needed to balance net global salt flux into ocean at surface', &
Expand Down
82 changes: 74 additions & 8 deletions src/parameterizations/vertical/MOM_diabatic_aux.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -67,7 +67,10 @@ module MOM_diabatic_aux
!! e-folding depth of incoming shortwave radiation.
type(external_field) :: sbc_chl !< A handle used in time interpolation of
!! chlorophyll read from a file.
logical :: chl_from_file !< If true, chl_a is read from a file.
logical :: chl_from_file !< If true, chl_a is read from a file.
logical :: do_brine_plume !< If true, insert salt flux below the surface according to
!! a parameterization by \cite Nguyen2009.
integer :: brine_plume_n !< The exponent in the brine plume parameterization.

type(time_type), pointer :: Time => NULL() !< A pointer to the ocean model's clock.
type(diag_ctrl), pointer :: diag !< Structure used to regulate timing of diagnostic output
Expand Down Expand Up @@ -1034,7 +1037,7 @@ end subroutine diagnoseMLDbyEnergy
subroutine applyBoundaryFluxesInOut(CS, G, GV, US, dt, fluxes, optics, nsw, h, tv, &
aggregate_FW_forcing, evap_CFL_limit, &
minimum_forcing_depth, cTKE, dSV_dT, dSV_dS, &
SkinBuoyFlux )
SkinBuoyFlux, MLD)
type(diabatic_aux_CS), pointer :: CS !< Control structure for diabatic_aux
type(ocean_grid_type), intent(in) :: G !< Grid structure
type(verticalGrid_type), intent(in) :: GV !< ocean vertical grid structure
Expand Down Expand Up @@ -1064,6 +1067,7 @@ subroutine applyBoundaryFluxesInOut(CS, G, GV, US, dt, fluxes, optics, nsw, h, t
!! salinity [R-1 S-1 ~> m3 kg-1 ppt-1].
real, dimension(SZI_(G),SZJ_(G)), &
optional, intent(out) :: SkinBuoyFlux !< Buoyancy flux at surface [Z2 T-3 ~> m2 s-3].
real, pointer, dimension(:,:), optional :: MLD!< Mixed layer depth for brine plumes [Z ~> m]

! Local variables
integer, parameter :: maxGroundings = 5
Expand Down Expand Up @@ -1102,7 +1106,8 @@ subroutine applyBoundaryFluxesInOut(CS, G, GV, US, dt, fluxes, optics, nsw, h, t
netheat_rate, & ! netheat but for dt=1 [C H T-1 ~> degC m s-1 or degC kg m-2 s-1]
netsalt_rate, & ! netsalt but for dt=1 (e.g. returns a rate)
! [S H T-1 ~> ppt m s-1 or ppt kg m-2 s-1]
netMassInOut_rate! netmassinout but for dt=1 [H T-1 ~> m s-1 or kg m-2 s-1]
netMassInOut_rate, & ! netmassinout but for dt=1 [H T-1 ~> m s-1 or kg m-2 s-1]
mixing_depth ! Mixed layer depth [Z -> m]
real, dimension(SZI_(G), SZK_(GV)) :: &
h2d, & ! A 2-d copy of the thicknesses [H ~> m or kg m-2]
T2d, & ! A 2-d copy of the layer temperatures [C ~> degC]
Expand Down Expand Up @@ -1132,13 +1137,19 @@ subroutine applyBoundaryFluxesInOut(CS, G, GV, US, dt, fluxes, optics, nsw, h, t
! and rejected brine are initially applied in vanishingly thin layers at the
! top of the layer before being mixed throughout the layer.
logical :: calculate_buoyancy ! If true, calculate the surface buoyancy flux.
real, dimension(SZI_(G)) :: dK ! Depth [Z ~> m].
real, dimension(SZI_(G)) :: A_brine ! Constant [Z-(n+1) ~> m-(n+1)].
real :: fraction_left_brine ! Sum for keeping track of the fraction of brine so far (in depth)
real :: plume_fraction ! Sum for keeping track of the fraction of brine so far (in depth)
real :: plume_flux ! Brine flux to move downwards [S H ~> ppt m or ppt kg m-2]
integer, dimension(2) :: EOSdom ! The i-computational domain for the equation of state
integer :: i, j, is, ie, js, je, k, nz, nb
character(len=45) :: mesg

is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke

Idt = 1.0 / dt
plume_flux = 0.0

calculate_energetics = (present(cTKE) .and. present(dSV_dT) .and. present(dSV_dS))
calculate_buoyancy = present(SkinBuoyFlux)
Expand All @@ -1158,6 +1169,17 @@ subroutine applyBoundaryFluxesInOut(CS, G, GV, US, dt, fluxes, optics, nsw, h, t
GoRho = US%L_to_Z**2*GV%g_Earth / GV%Rho0
endif

if (CS%do_brine_plume .and. .not. associated(MLD)) then
call MOM_error(FATAL, "MOM_diabatic_aux.F90, applyBoundaryFluxesInOut(): "//&
"Brine plume parameterization requires a mixed-layer depth,\n"//&
"currently coming from the energetic PBL scheme.")
endif
if (CS%do_brine_plume .and. .not. associated(fluxes%salt_left_behind)) then
call MOM_error(FATAL, "MOM_diabatic_aux.F90, applyBoundaryFluxesInOut(): "//&
"Brine plume parameterization requires DO_BRINE_PLUME\n"//&
"to be turned on in SIS2 as well as MOM6.")
endif

! H_limit_fluxes is used by extractFluxes1d to scale down fluxes if the total
! depth of the ocean is vanishing. It does not (yet) handle a value of zero.
! To accommodate vanishing upper layers, we need to allow for an instantaneous
Expand All @@ -1173,17 +1195,19 @@ subroutine applyBoundaryFluxesInOut(CS, G, GV, US, dt, fluxes, optics, nsw, h, t
!$OMP H_limit_fluxes,numberOfGroundings,iGround,jGround,&
!$OMP nonPenSW,hGrounding,CS,Idt,aggregate_FW_forcing, &
!$OMP minimum_forcing_depth,evap_CFL_limit,dt,EOSdom, &
!$OMP calculate_buoyancy,netPen_rate,SkinBuoyFlux,GoRho, &
!$OMP calculate_buoyancy,netPen_rate,SkinBuoyFlux,GoRho,&
!$OMP calculate_energetics,dSV_dT,dSV_dS,cTKE,g_Hconv2, &
!$OMP EnthalpyConst) &
!$OMP EnthalpyConst,MLD) &
!$OMP private(opacityBand,h2d,T2d,netMassInOut,netMassOut, &
!$OMP netHeat,netSalt,Pen_SW_bnd,fractionOfForcing, &
!$OMP IforcingDepthScale, &
!$OMP dThickness,dTemp,dSalt,hOld,Ithickness, &
!$OMP netMassIn,pres,d_pres,p_lay,dSV_dT_2d, &
!$OMP netmassinout_rate,netheat_rate,netsalt_rate, &
!$OMP drhodt,drhods,pen_sw_bnd_rate, &
!$OMP pen_TKE_2d,Temp_in,Salin_in,RivermixConst) &
!$OMP pen_TKE_2d,Temp_in,Salin_in,RivermixConst, &
!$OMP mixing_depth,A_brine,fraction_left_brine, &
!$OMP plume_flux,plume_fraction,dK) &
!$OMP firstprivate(SurfPressure)
do j=js,je
! Work in vertical slices for efficiency
Expand Down Expand Up @@ -1300,6 +1324,14 @@ subroutine applyBoundaryFluxesInOut(CS, G, GV, US, dt, fluxes, optics, nsw, h, t
! ocean (and corresponding outward heat content), and ignoring penetrative SW.
! B/ update mass, salt, temp from mass leaving ocean.
! C/ update temp due to penetrative SW
if (CS%do_brine_plume) then
do i=is,ie
mixing_depth(i) = max(MLD(i,j) - minimum_forcing_depth * GV%H_to_Z, minimum_forcing_depth * GV%H_to_Z)
mixing_depth(i) = min(mixing_depth(i), max(sum(h(i,j,:)), GV%angstrom_h) * GV%H_to_Z)
A_brine(i) = (CS%brine_plume_n + 1) / (mixing_depth(i) ** (CS%brine_plume_n + 1))
enddo
endif

do i=is,ie
if (G%mask2dT(i,j) > 0.) then

Expand Down Expand Up @@ -1372,8 +1404,8 @@ subroutine applyBoundaryFluxesInOut(CS, G, GV, US, dt, fluxes, optics, nsw, h, t
enddo ! k=1,1

! B/ Update mass, salt, temp from mass leaving ocean and other fluxes of heat and salt.
fraction_left_brine = 1.0
do k=1,nz

! Place forcing into this layer if this layer has nontrivial thickness.
! For layers thin relative to 1/IforcingDepthScale, then distribute
! forcing into deeper layers.
Expand All @@ -1388,6 +1420,33 @@ subroutine applyBoundaryFluxesInOut(CS, G, GV, US, dt, fluxes, optics, nsw, h, t
fractionOfForcing = -evap_CFL_limit*h2d(i,k)/netMassOut(i)
endif

if (CS%do_brine_plume .and. associated(fluxes%salt_left_behind)) then
if (fluxes%salt_left_behind(i,j) > 0 .and. fraction_left_brine > 0.0) then
! Place forcing into this layer by depth for brine plume parameterization.
if (k == 1) then
dK(i) = 0.5 * h(i,j,k) * GV%H_to_Z ! Depth of center of layer K
plume_flux = - (1000.0*US%ppt_to_S * fluxes%salt_left_behind(i,j)) * GV%RZ_to_H
plume_fraction = 1.0
else
dK(i) = dK(i) + 0.5 * ( h(i,j,k) + h(i,j,k-1) ) * GV%H_to_Z ! Depth of center of layer K
plume_flux = 0.0
endif
if (dK(i) <= mixing_depth(i) .and. fraction_left_brine > 0.0) then
plume_fraction = min(fraction_left_brine, A_brine(i) * dK(i) ** CS%brine_plume_n &
* h(i,j,k) * GV%H_to_Z)
else
IforcingDepthScale = 1. / max(GV%H_subroundoff, minimum_forcing_depth - netMassOut(i) )
! plume_fraction = fraction_left_brine, unless h2d is less than IforcingDepthScale.
plume_fraction = min(fraction_left_brine, h2d(i,k)*IforcingDepthScale)
endif
fraction_left_brine = fraction_left_brine - plume_fraction
plume_flux = plume_flux + plume_fraction * (1000.0*US%ppt_to_S * fluxes%salt_left_behind(i,j)) &
* GV%RZ_to_H
else
plume_flux = 0.0
endif
endif

! Change in state due to forcing

dThickness = max( fractionOfForcing*netMassOut(i), -h2d(i,k) )
Expand Down Expand Up @@ -1432,7 +1491,7 @@ subroutine applyBoundaryFluxesInOut(CS, G, GV, US, dt, fluxes, optics, nsw, h, t
endif
Ithickness = 1.0/h2d(i,k) ! Inverse of new thickness
T2d(i,k) = (hOld*T2d(i,k) + dTemp)*Ithickness
tv%S(i,j,k) = (hOld*tv%S(i,j,k) + dSalt)*Ithickness
tv%S(i,j,k) = (hOld*tv%S(i,j,k) + dSalt + plume_flux)*Ithickness
elseif (h2d(i,k) < 0.0) then ! h2d==0 is a special limit that needs no extra handling
call forcing_SinglePointPrint(fluxes,G,i,j,'applyBoundaryFluxesInOut (h<0)')
write(0,*) 'applyBoundaryFluxesInOut(): lon,lat=',G%geoLonT(i,j),G%geoLatT(i,j)
Expand Down Expand Up @@ -1702,6 +1761,13 @@ subroutine diabatic_aux_init(Time, G, GV, US, param_file, diag, CS, useALEalgori
CS%use_calving_heat_content = .false.
endif

call get_param(param_file, mdl, "DO_BRINE_PLUME", CS%do_brine_plume, &
"If true, use a brine plume parameterization from "//&
"Nguyen et al., 2009.", default=.false.)
call get_param(param_file, mdl, "BRINE_PLUME_EXPONENT", CS%brine_plume_n, &
"If using the brine plume parameterization, set the integer exponent.", &
default=5, do_not_log=.not.CS%do_brine_plume)

if (useALEalgorithm) then
CS%id_createdH = register_diag_field('ocean_model',"created_H",diag%axesT1, &
Time, "The volume flux added to stop the ocean from drying out and becoming negative in depth", &
Expand Down
Loading

0 comments on commit 147ddf1

Please sign in to comment.