Skip to content

Commit

Permalink
Merge pull request #548 from breichl/user/bgr/Fix_BuoyFlux_merged
Browse files Browse the repository at this point in the history 
User/bgr/fix buoy flux merged
  • Loading branch information
@adcroft
adcroft authored Jul 7, 2017
2 parents da325d6 + 8a5704d commit ea1a338
Show file tree
Hide file tree
Showing 3 changed files with 187 additions and 20 deletions.
87 changes: 85 additions & 2 deletions src/core/MOM_forcing_type.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -273,7 +273,8 @@ module MOM_forcing_type
subroutine extractFluxes1d(G, GV, fluxes, optics, nsw, j, dt, &
DepthBeforeScalingFluxes, useRiverHeatContent, useCalvingHeatContent, &
h, T, netMassInOut, netMassOut, net_heat, net_salt, pen_SW_bnd, tv, &
aggregate_FW_forcing, nonpenSW, skip_diags)
aggregate_FW_forcing, nonpenSW, netmassInOut_rate,net_Heat_Rate, &
net_salt_rate, pen_sw_bnd_Rate, skip_diags)

type(ocean_grid_type), intent(in) :: G !< ocean grid structure
type(verticalGrid_type), intent(in) :: GV !< ocean vertical grid structure
Expand Down Expand Up @@ -317,12 +318,17 @@ subroutine extractFluxes1d(G, GV, fluxes, optics, nsw, j, dt,
real, dimension(SZI_(G)), optional, intent(out) :: nonpenSW !< non-downwelling SW; use in net_heat.
!! Sum over SW bands when diagnosing nonpenSW.
!! Units are (K * H).
real, dimension(SZI_(G)), optional, intent(out) :: net_Heat_rate !< Optional outputs of contributions to surface
real, dimension(SZI_(G)), optional, intent(out) :: net_salt_rate !< buoyancy flux which do not include dt
real, dimension(SZI_(G)), optional, intent(out) :: netmassInOut_rate !< and therefore are used to compute the rate.
real, dimension(:,:), optional, intent(out) :: pen_sw_bnd_rate !< Perhaps just a temporary fix.
logical, optional, intent(in) :: skip_diags !< If present and true, skip calculating
!! diagnostics

! local
real :: htot(SZI_(G)) ! total ocean depth (m for Bouss or kg/m^2 for non-Bouss)
real :: Pen_sw_tot(SZI_(G)) ! sum across all bands of Pen_SW (K * H)
real :: pen_sw_tot_rate(SZI_(G)) ! Similar but sum but as a rate (no dt in calculation)
real :: Ih_limit ! inverse depth at which surface fluxes start to be limited (1/H)
real :: scale ! scale scales away fluxes if depth < DepthBeforeScalingFluxes
real :: J_m2_to_H ! converts J/m^2 to H units (m for Bouss and kg/m^2 for non-Bouss)
Expand All @@ -331,6 +337,23 @@ subroutine extractFluxes1d(G, GV, fluxes, optics, nsw, j, dt,
logical :: calculate_diags ! Indicate to calculate/update diagnostic arrays
character(len=200) :: mesg
integer :: is, ie, nz, i, k, n

logical :: do_NHR, do_NSR, do_NMIOR, do_PSWBR

!BGR-Jul 5,2017{
! Initializes/sets logicals if 'rates' are requested
! These factors are required for legacy reasons
! and therefore computed only when optional outputs are requested
do_NHR = .false.
do_NSR = .false.
do_NMIOR = .false.
do_PSWBR = .false.
if (present(net_heat_rate)) do_NHR = .true.
if (present(net_salt_rate)) do_NSR = .true.
if (present(netmassinout_rate)) do_NMIOR = .true.
if (present(pen_sw_bnd_rate)) do_PSWBR = .true.
!}BGR

Ih_limit = 1.0 / DepthBeforeScalingFluxes
Irho0 = 1.0 / GV%Rho0
I_Cp = 1.0 / fluxes%C_p
Expand Down Expand Up @@ -390,6 +413,21 @@ subroutine extractFluxes1d(G, GV, fluxes, optics, nsw, j, dt,
Pen_SW_bnd(1,i) = 0.0
endif

!BGR-Jul 5, 2017{
!Repeats above code w/ dt=1. for legacy reason
if (do_PSWBR) then
pen_sw_tot_rate(i) = 0.0
if (nsw >= 1) then
do n=1,nsw
Pen_SW_bnd_rate(n,i) = J_m2_to_H*scale * max(0.0, optics%sw_pen_band(n,i,j))
pen_sw_tot_rate(i) = pen_sw_tot_rate(i) + pen_sw_bnd_rate(n,i)
enddo
else
pen_sw_bnd_rate(1,i) = 0.0
endif
endif
!}BGR

! net volume/mass of liquid and solid passing through surface boundary fluxes
netMassInOut(i) = dt * (scale * ((((( fluxes%lprec(i,j) &
+ fluxes%fprec(i,j) ) &
Expand All @@ -398,12 +436,29 @@ subroutine extractFluxes1d(G, GV, fluxes, optics, nsw, j, dt,
+ fluxes%vprec(i,j) ) &
+ fluxes%frunoff(i,j) ) )

!BGR-Jul 5, 2017{
!Repeats above code w/ dt=1. for legacy reason
if (do_NMIOr) then
netMassInOut_rate(i) = (scale * ((((( fluxes%lprec(i,j) &
+ fluxes%fprec(i,j) ) &
+ fluxes%evap(i,j) ) &
+ fluxes%lrunoff(i,j) ) &
+ fluxes%vprec(i,j) ) &
+ fluxes%frunoff(i,j) ) )
endif
!}BGR

! smg:
! for non-Bouss, we add/remove salt mass to total ocean mass. to conserve
! total salt mass ocean+ice, the sea ice model must lose mass when
! salt mass is added to the ocean, which may still need to be coded.
if (.not.GV%Boussinesq .and. ASSOCIATED(fluxes%salt_flux)) then
netMassInOut(i) = netMassInOut(i) + (dt * GV%kg_m2_to_H) * (scale * fluxes%salt_flux(i,j))

!BGR-Jul 5, 2017{
!Repeats above code w/ dt=1. for legacy reason
if (do_NMIOr) netMassInOut_rate(i) = netMassInOut_rate(i) + (GV%kg_m2_to_H) * (scale * fluxes%salt_flux(i,j))
!}BGR
endif

! net volume/mass of water leaving the ocean.
Expand Down Expand Up @@ -433,16 +488,25 @@ subroutine extractFluxes1d(G, GV, fluxes, optics, nsw, j, dt,

! convert to H units (Bouss=meter or non-Bouss=kg/m^2)
netMassInOut(i) = GV%kg_m2_to_H * netMassInOut(i)
!BGR-Jul 5, 2017{
!Repeats above code w/ dt=1. for legacy reason
if (do_NMIOr) netMassInOut_rate(i) = GV%kg_m2_to_H * netMassInOut_rate(i)
!}BGR
netMassOut(i) = GV%kg_m2_to_H * netMassOut(i)

! surface heat fluxes from radiation and turbulent fluxes (K * H)
! (H=m for Bouss, H=kg/m2 for non-Bouss)
net_heat(i) = scale * dt * J_m2_to_H * &
( fluxes%sw(i,j) + ((fluxes%lw(i,j) + fluxes%latent(i,j)) + fluxes%sens(i,j)) )

!BGR-Jul 5, 2017{
!Repeats above code w/ dt=1. for legacy reason
if (do_NHR) net_heat_rate(i) = scale * J_m2_to_H * &
( fluxes%sw(i,j) + ((fluxes%lw(i,j) + fluxes%latent(i,j)) + fluxes%sens(i,j)) )
!}BGR
! Add heat flux from surface damping (restoring) (K * H) or flux adjustments.
if (ASSOCIATED(fluxes%heat_added)) then
net_heat(i) = net_heat(i) + (scale * (dt * J_m2_to_H)) * fluxes%heat_added(i,j)
if (do_NHR) net_heat_rate(i) = net_heat_rate(i) + (scale * (J_m2_to_H)) * fluxes%heat_added(i,j)
endif

! Add explicit heat flux for runoff (which is part of the ice-ocean boundary
Expand All @@ -451,6 +515,11 @@ subroutine extractFluxes1d(G, GV, fluxes, optics, nsw, j, dt,
! remove lrunoff*SST here, to counteract its addition elsewhere
net_heat(i) = (net_heat(i) + (scale*(dt*J_m2_to_H)) * fluxes%heat_content_lrunoff(i,j)) - &
(GV%kg_m2_to_H * (scale * dt)) * fluxes%lrunoff(i,j) * T(i,1)
!BGR-Jul 5, 2017{
!Intentionally neglect the following contribution to rate for legacy reasons.
!if (do_NHR) net_heat_rate(i) = (net_heat_rate(i) + (scale*(J_m2_to_H)) * fluxes%heat_content_lrunoff(i,j)) - &
! (GV%kg_m2_to_H * (scale)) * fluxes%lrunoff(i,j) * T(i,1)
!}BGR
if (calculate_diags .and. ASSOCIATED(tv%TempxPmE)) then
tv%TempxPmE(i,j) = tv%TempxPmE(i,j) + (scale * dt) * &
(I_Cp*fluxes%heat_content_lrunoff(i,j) - fluxes%lrunoff(i,j)*T(i,1))
Expand All @@ -463,6 +532,11 @@ subroutine extractFluxes1d(G, GV, fluxes, optics, nsw, j, dt,
! remove frunoff*SST here, to counteract its addition elsewhere
net_heat(i) = net_heat(i) + (scale*(dt*J_m2_to_H)) * fluxes%heat_content_frunoff(i,j) - &
(GV%kg_m2_to_H * (scale * dt)) * fluxes%frunoff(i,j) * T(i,1)
!BGR-Jul 5, 2017{
!Intentionally neglect the following contribution to rate for legacy reasons.
! if (do_NHR) net_heat_rate(i) = net_heat_rate(i) + (scale*(J_m2_to_H)) * fluxes%heat_content_frunoff(i,j) - &
! (GV%kg_m2_to_H * (scale)) * fluxes%frunoff(i,j) * T(i,1)
!}BGR
if (calculate_diags .and. ASSOCIATED(tv%TempxPmE)) then
tv%TempxPmE(i,j) = tv%TempxPmE(i,j) + (scale * dt) * &
(I_Cp*fluxes%heat_content_frunoff(i,j) - fluxes%frunoff(i,j)*T(i,1))
Expand Down Expand Up @@ -500,6 +574,10 @@ subroutine extractFluxes1d(G, GV, fluxes, optics, nsw, j, dt,
! remove penetrative portion of the SW that is NOT absorbed within a
! tiny layer at the top of the ocean.
net_heat(i) = net_heat(i) - Pen_SW_tot(i)
!BGR-Jul 5, 2017{
!Repeat above code for 'rate' term
if (do_NHR) net_heat_rate(i) = net_heat_rate(i) - Pen_SW_tot_rate(i)
!}BGR

! diagnose non-downwelling SW
if (present(nonPenSW)) then
Expand All @@ -508,11 +586,16 @@ subroutine extractFluxes1d(G, GV, fluxes, optics, nsw, j, dt,

! Salt fluxes
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 * fluxes%salt_flux(i,j))) * GV%kg_m2_to_H
!BGR-Jul 5, 2017{
!Repeat above code for 'rate' term
if (do_NSR) Net_salt_rate(i) = (scale * 1. * (1000.0 * fluxes%salt_flux(i,j))) * GV%kg_m2_to_H
!}BGR
endif

! Diagnostics follow...
Expand Down
101 changes: 93 additions & 8 deletions src/parameterizations/vertical/MOM_diabatic_aux.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -71,14 +71,14 @@ module MOM_diabatic_aux
use MOM_diag_mediator, only : post_data, register_diag_field, safe_alloc_ptr
use MOM_diag_mediator, only : diag_ctrl, time_type
use MOM_EOS, only : calculate_density, calculate_TFreeze
use MOM_EOS, only : calculate_specific_vol_derivs
use MOM_EOS, only : calculate_specific_vol_derivs, calculate_density_derivs
use MOM_error_handler, only : MOM_error, FATAL, WARNING, callTree_showQuery
use MOM_error_handler, only : callTree_enter, callTree_leave, callTree_waypoint
use MOM_file_parser, only : get_param, log_version, param_file_type
use MOM_forcing_type, only : forcing, extractFluxes1d, forcing_SinglePointPrint
use MOM_grid, only : ocean_grid_type
use MOM_io, only : vardesc
use MOM_shortwave_abs, only : absorbRemainingSW, optics_type
use MOM_shortwave_abs, only : absorbRemainingSW, optics_type, sumSWoverBands
use MOM_variables, only : thermo_var_ptrs, vertvisc_type! , accel_diag_ptrs
use MOM_verticalGrid, only : verticalGrid_type

Expand Down Expand Up @@ -809,7 +809,8 @@ end subroutine diagnoseMLDbyDensityDifference
!! and calculate the TKE implications of this heating.
subroutine applyBoundaryFluxesInOut(CS, G, GV, dt, fluxes, optics, h, tv, &
aggregate_FW_forcing, evap_CFL_limit, &
minimum_forcing_depth, cTKE, dSV_dT, dSV_dS)
minimum_forcing_depth, cTKE, dSV_dT, dSV_dS, &
SkinBuoyFlux )
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 All @@ -830,6 +831,8 @@ subroutine applyBoundaryFluxesInOut(CS, G, GV, dt, fluxes, optics, h, tv, &
real, dimension(SZI_(G),SZJ_(G),SZK_(G)), optional, intent(out) :: dSV_dT
!> Partial derivative of specific a volume with potential salinity, in m3 kg-1 / (g kg-1).
real, dimension(SZI_(G),SZJ_(G),SZK_(G)), optional, intent(out) :: dSV_dS
!> Buoyancy flux at surface in m2 s-3
real, dimension(SZI_(G),SZJ_(G)), optional, intent(out) :: SkinBuoyFlux

! Local variables
integer, parameter :: maxGroundings = 5
Expand All @@ -848,16 +851,27 @@ subroutine applyBoundaryFluxesInOut(CS, G, GV, dt, fluxes, optics, h, tv, &
netHeat, & ! heat (degC * H) via surface fluxes, excluding
! Pen_SW_bnd and netMassOut
netSalt, & ! surface salt flux ( g(salt)/m2 for non-Bouss and ppt*H for Bouss )
nonpenSW ! non-downwelling SW, which is absorbed at ocean surface
nonpenSW, & ! non-downwelling SW, which is absorbed at ocean surface
SurfPressure, & ! Surface pressure (approximated as 0.0)
dRhodT, & ! change in density per change in temperature
dRhodS, & ! change in density per change in salinity
netheat_rate, & ! netheat but for dt=1 (e.g. returns a rate)
netsalt_rate, & ! netsalt but for dt=1 (e.g. returns a rate)
netMassInOut_rate! netmassinout but for dt=1 (e.g. returns a rate)
real, dimension(SZI_(G), SZK_(G)) :: h2d, T2d
real, dimension(SZI_(G), SZK_(G)) :: pen_TKE_2d, dSV_dT_2d
real, dimension(max(optics%nbands,1),SZI_(G)) :: Pen_SW_bnd
real, dimension(SZI_(G),SZK_(G)+1) :: netPen
real, dimension(max(optics%nbands,1),SZI_(G)) :: Pen_SW_bnd, Pen_SW_bnd_rate
!^ _rate is w/ dt=1
real, dimension(max(optics%nbands,1),SZI_(G),SZK_(G)) :: opacityBand
real :: hGrounding(maxGroundings)
real :: Temp_in, Salin_in
real :: I_G_Earth, g_Hconv2
real :: GoRho
logical :: calculate_energetics
logical :: calculate_buoyancy
integer :: i, j, is, ie, js, je, k, nz, n, nsw
integer :: start, npts
character(len=45) :: mesg

is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = G%ke
Expand All @@ -870,10 +884,18 @@ subroutine applyBoundaryFluxesInOut(CS, G, GV, dt, fluxes, optics, h, tv, &
Idt = 1.0/dt

calculate_energetics = (present(cTKE) .and. present(dSV_dT) .and. present(dSV_dS))
calculate_buoyancy = present(SkinBuoyFlux)
if (calculate_buoyancy) SkinBuoyFlux(:,:) = 0.0
I_G_Earth = 1.0 / GV%g_Earth
g_Hconv2 = GV%g_Earth * GV%H_to_kg_m2**2

if (present(cTKE)) cTKE(:,:,:) = 0.0
if (calculate_buoyancy) then
SurfPressure(:) = 0.0
GoRho = GV%g_Earth / GV%Rho0
start = 1 + G%isc - G%isd
npts = 1 + G%iec - G%isc
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.
Expand Down Expand Up @@ -947,11 +969,47 @@ subroutine applyBoundaryFluxesInOut(CS, G, GV, dt, fluxes, optics, h, tv, &
! enters to the ocean and participates in pentrative SW heating.
! nonpenSW = non-downwelling SW flux, which is absorbed in ocean surface
! (in tandem w/ LW,SENS,LAT); saved only for diagnostic purposes.
call extractFluxes1d(G, GV, fluxes, optics, nsw, j, dt, &

!----------------------------------------------------------------------------------------
!BGR-June 26, 2017{
!Temporary action to preserve answers while fixing a bug.
! To fix a bug in a diagnostic calculation, applyboundaryfluxesinout now returns
! the surface buoyancy flux. Previously, extractbuoyancyflux2d was called, meaning
! a second call to extractfluxes1d (causing the diagnostic net_heat to be incorrect).
! Note that this call to extract buoyancyflux2d was AFTER applyboundaryfluxesinout,
! which means it used the T/S fields after this routine. Therefore, the surface
! buoyancy flux is computed here at the very end of this routine for legacy reasons.
! A few specific notes follow:
! 1) The old method did not included river/calving contributions to heat flux. This
! is kept consistent here via commenting code in the present extractFluxes1d <_rate>
! outputs, but we may reconsider this approach.
! 2) The old method computed the buoyancy flux rate directly (by setting dt=1), instead
! of computing the integrated value (and dividing by dt). Hence the required
! additional outputs from extractFluxes1d.
! *** This is because: A*dt/dt =/= A due to round off.
! 3) The old method computed buoyancy flux after this routine, meaning the returned
! surface fluxes (from extractfluxes1d) must be recorded for use later in the code.
! We could (and maybe should) move that loop up to before the surface fluxes are
! applied, but this will change answers.
! For all these reasons we compute additional values of <_rate> which are preserved
! for the buoyancy flux calculation and reproduce the old answers.
! In the future this needs more detailed investigation to make sure everything is
! consistent and correct. These details shouldnt significantly effect climate,
! but do change answers.
!-----------------------------------------------------------------------------------------
if (calculate_buoyancy) then
call extractFluxes1d(G, GV, fluxes, optics, nsw, j, dt, &
H_limit_fluxes, CS%use_river_heat_content, CS%use_calving_heat_content, &
h2d, T2d, netMassInOut, netMassOut, netHeat, netSalt, &
Pen_SW_bnd, tv, aggregate_FW_forcing, nonpenSW)

Pen_SW_bnd, tv, aggregate_FW_forcing, nonpenSW=nonpenSW, &
net_Heat_rate=netheat_rate,net_salt_rate=netsalt_rate, &
netmassinout_rate=netmassinout_rate,pen_sw_bnd_rate=pen_sw_bnd_rate)
else
call extractFluxes1d(G, GV, fluxes, optics, nsw, j, dt, &
H_limit_fluxes, CS%use_river_heat_content, CS%use_calving_heat_content, &
h2d, T2d, netMassInOut, netMassOut, netHeat, netSalt, &
Pen_SW_bnd, tv, aggregate_FW_forcing, nonpenSW=nonpenSW)
endif
! ea is for passive tracers
do i=is,ie
! ea(i,j,1) = netMassInOut(i)
Expand Down Expand Up @@ -1211,6 +1269,33 @@ subroutine applyBoundaryFluxesInOut(CS, G, GV, dt, fluxes, optics, h, tv, &
enddo
endif

! BGR: Get buoyancy flux to return for ePBL
! We want the rate, so we use the rate values returned from extractfluxes1d.
! Note that the *dt values could be divided by dt here, but
! 1) Answers will change due to round-off
! 2) Be sure to save their values BEFORE fluxes are used.
if (Calculate_Buoyancy) then
drhodt(:) = 0.0
drhods(:) = 0.0
netPen(:,:) = 0.0
! Sum over bands and attenuate as a function of depth
! netPen is the netSW as a function of depth
call sumSWoverBands(G, GV, h2d(:,:), optics%opacity_band(:,:,j,:), nsw, j, dt, &
H_limit_fluxes, .true., pen_SW_bnd_rate, netPen)
! Density derivatives
call calculate_density_derivs(T2d(:,1), tv%S(:,j,1), SurfPressure, &
dRhodT, dRhodS, start, npts, tv%eqn_of_state)
! 1. Adjust netSalt to reflect dilution effect of FW flux
! 2. Add in the SW heating for purposes of calculating the net
! surface buoyancy flux affecting the top layer.
! 3. Convert to a buoyancy flux, excluding penetrating SW heating
! BGR-Jul 5, 2017: The contribution of SW heating here needs investigated for ePBL.
SkinBuoyFlux(G%isc:G%iec,j) = - GoRho * ( dRhodS(G%isc:G%iec) * (netSalt_rate(G%isc:G%iec) &
- tv%S(G%isc:G%iec,j,1) * netMassInOut_rate(G%isc:G%iec)* GV%H_to_m )&
+ dRhodT(G%isc:G%iec) * ( netHeat_rate(G%isc:G%iec) + &
netPen(G%isc:G%iec,1))) * GV%H_to_m ! m^2/s^3
endif

enddo ! j-loop finish

! Post the diagnostics
Expand Down
Loading

0 comments on commit ea1a338

Please sign in to comment.