Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

User/bgr/fix buoy flux merged #548

Merged
merged 6 commits into from
Jul 7, 2017
Merged
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
87 changes: 85 additions & 2 deletions src/core/MOM_forcing_type.F90
Original file line number Diff line number Diff line change
@@ -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
@@ -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)
@@ -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
@@ -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) ) &
@@ -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.
@@ -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
@@ -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))
@@ -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))
@@ -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
@@ -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...
101 changes: 93 additions & 8 deletions src/parameterizations/vertical/MOM_diabatic_aux.F90
Original file line number Diff line number Diff line change
@@ -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

@@ -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
@@ -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
@@ -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
@@ -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.
@@ -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)
@@ -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
Loading