Skip to content

Commit

Permalink
+Rescaled _rate variables from extractFluxes1d
Browse files Browse the repository at this point in the history 
  Rescaled the units of optional _rate variables returned by extractFluxes1d and
simplified the calculations using these variables in applyBoundaryFluxesInOut.
All answers are bitwise identical, but the units of arguments to a public type
have changed.
  • Loading branch information
@Hallberg-NOAA
Hallberg-NOAA committed Oct 5, 2019
1 parent 76172d4 commit 965c2c2
Show file tree
Hide file tree
Showing 2 changed files with 39 additions and 30 deletions.
35 changes: 19 additions & 16 deletions src/core/MOM_forcing_type.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -339,7 +339,7 @@ module MOM_forcing_type
subroutine extractFluxes1d(G, GV, US, fluxes, optics, nsw, j, dt, &
FluxRescaleDepth, useRiverHeatContent, useCalvingHeatContent, &
h, T, netMassInOut, netMassOut, net_heat, net_salt, pen_SW_bnd, tv, &
aggregate_FW, nonpenSW, netmassInOut_rate,net_Heat_Rate, &
aggregate_FW, 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
Expand Down Expand Up @@ -392,22 +392,23 @@ subroutine extractFluxes1d(G, GV, US, fluxes, optics, nsw, j, dt,
!! Summed over SW bands when diagnosing nonpenSW.
real, dimension(SZI_(G)), &
optional, intent(out) :: net_Heat_rate !< Rate of net surface heating
!! [degC H s-1 ~> degC m s-1 or degC kg m-2 s-1].
!! [degC H T-1 ~> degC m s-1 or degC kg m-2 s-1].
real, dimension(SZI_(G)), &
optional, intent(out) :: net_salt_rate !< Surface salt flux into the ocean
!! [ppt H s-1 ~> ppt m s-1 or ppt kg m-2 s-1].
!! [ppt H T-1 ~> ppt m s-1 or ppt kg m-2 s-1].
real, dimension(SZI_(G)), &
optional, intent(out) :: netmassInOut_rate !< Rate of net mass flux into the ocean
!! [H s-1 ~> m s-1 or kg m-2 s-1].
!! [H T-1 ~> m s-1 or kg m-2 s-1].
real, dimension(max(1,nsw),G%isd:G%ied), &
optional, intent(out) :: pen_sw_bnd_rate !< Rate of penetrative shortwave heating
!! [degC H s-1 ~> degC m s-1 or degC kg m-2 s-1].
!! [degC H T-1 ~> degC m s-1 or degC kg m-2 s-1].
logical, optional, intent(in) :: skip_diags !< If present and true, skip calculating diagnostics

! local
real :: htot(SZI_(G)) ! total ocean depth [H ~> m or kg m-2]
real :: Pen_sw_tot(SZI_(G)) ! sum across all bands of Pen_SW [degC H ~> degC m or degC kg m-2].
real :: pen_sw_tot_rate(SZI_(G)) ! Similar but sum but as a rate (no dt in calculation)
real :: pen_sw_tot_rate(SZI_(G)) ! Summed rate of shortwave heating across bands
! [degC H T-1 ~> degC m s-1 or degC kg m-2 s-1]
real :: Ih_limit ! inverse depth at which surface fluxes start to be limited [H-1 ~> m-1 or m2 kg-1]
real :: scale ! scale scales away fluxes if depth < FluxRescaleDepth
real :: J_m2_to_H ! converts J/m^2 to H units (m for Bouss and kg/m^2 for non-Bouss)
Expand Down Expand Up @@ -503,7 +504,7 @@ subroutine extractFluxes1d(G, GV, US, fluxes, optics, nsw, j, dt,
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, Pen_SW_bnd_rate(n,i))
Pen_SW_bnd_rate(n,i) = J_m2_to_H*US%T_to_s*scale * max(0.0, Pen_SW_bnd_rate(n,i))
pen_sw_tot_rate(i) = pen_sw_tot_rate(i) + pen_sw_bnd_rate(n,i)
enddo
else
Expand All @@ -521,8 +522,8 @@ subroutine extractFluxes1d(G, GV, US, fluxes, optics, nsw, j, dt,
+ fluxes%seaice_melt(i,j)) &
+ fluxes%frunoff(i,j) ))

if (do_NMIOr) then ! Repeat the above code w/ dt=1s for legacy reasons
netMassInOut_rate(i) = (scale * US%s_to_T* &
if (do_NMIOr) then ! Repeat the above code without multiplying by a timestep for legacy reasons
netMassInOut_rate(i) = (scale * &
(((((( fluxes%lprec(i,j) &
+ fluxes%fprec(i,j) ) &
+ fluxes%evap(i,j) ) &
Expand All @@ -540,7 +541,7 @@ subroutine extractFluxes1d(G, GV, US, fluxes, optics, nsw, j, dt,
if (.not.GV%Boussinesq .and. associated(fluxes%salt_flux)) then
netMassInOut(i) = netMassInOut(i) + dt * (scale * US%kg_m3_to_R*US%m_to_Z*fluxes%salt_flux(i,j))
if (do_NMIOr) netMassInOut_rate(i) = netMassInOut_rate(i) + &
(scale * US%kg_m3_to_R*US%m_to_Z*fluxes%salt_flux(i,j))
(scale * US%kg_m3_to_R*US%m_to_Z*US%T_to_s*fluxes%salt_flux(i,j))
endif

! net volume/mass of water leaving the ocean.
Expand Down Expand Up @@ -580,21 +581,21 @@ subroutine extractFluxes1d(G, GV, US, fluxes, optics, nsw, j, dt,
( fluxes%sw(i,j) + ((fluxes%lw(i,j) + fluxes%latent(i,j)) + fluxes%sens(i,j) + &
fluxes%seaice_melt_heat(i,j)) )
!Repeats above code w/ dt=1. for legacy reason
if (do_NHR) net_heat_rate(i) = scale * J_m2_to_H * &
if (do_NHR) net_heat_rate(i) = scale * US%T_to_s*J_m2_to_H * &
( fluxes%sw(i,j) + ((fluxes%lw(i,j) + fluxes%latent(i,j)) + fluxes%sens(i,j) + &
fluxes%seaice_melt_heat(i,j)))
else
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)) )
!Repeats above code w/ dt=1. for legacy reason
if (do_NHR) net_heat_rate(i) = scale * J_m2_to_H * &
if (do_NHR) net_heat_rate(i) = scale * US%T_to_s*J_m2_to_H * &
( fluxes%sw(i,j) + ((fluxes%lw(i,j) + fluxes%latent(i,j)) + fluxes%sens(i,j)) )
endif

! 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)
if (do_NHR) net_heat_rate(i) = net_heat_rate(i) + (scale * (US%T_to_s*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 @@ -605,7 +606,8 @@ subroutine extractFluxes1d(G, GV, US, fluxes, optics, nsw, j, dt,
(GV%RZ_to_H * (scale * dt_in_T)) * 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)) - &
!if (do_NHR) net_heat_rate(i) = (net_heat_rate(i) + (scale*(US%T_to_s*J_m2_to_H)) * &
! fluxes%heat_content_lrunoff(i,j)) - &
! (GV%RZ_to_H * (scale)) * US%s_to_T*fluxes%lrunoff(i,j) * T(i,1)
!}BGR
if (calculate_diags .and. associated(tv%TempxPmE)) then
Expand All @@ -622,7 +624,8 @@ subroutine extractFluxes1d(G, GV, US, fluxes, optics, nsw, j, dt,
(GV%RZ_to_H * (scale * dt_in_T)) * 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) - &
! if (do_NHR) net_heat_rate(i) = net_heat_rate(i) + (scale*(US%T_to_s*J_m2_to_H)) * &
! fluxes%heat_content_frunoff(i,j) - &
! (GV%RZ_to_H * (scale)) * US%s_to_T*fluxes%frunoff(i,j) * T(i,1)
!}BGR
if (calculate_diags .and. associated(tv%TempxPmE)) then
Expand Down Expand Up @@ -679,7 +682,7 @@ subroutine extractFluxes1d(G, GV, US, fluxes, optics, nsw, j, dt,
if (associated(fluxes%salt_flux)) then
Net_salt(i) = (scale * dt * (1000.0 * fluxes%salt_flux(i,j))) * GV%kg_m2_to_H
!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
if (do_NSR) Net_salt_rate(i) = (scale * 1. * (1000.0 * US%T_to_s*fluxes%salt_flux(i,j))) * GV%kg_m2_to_H
endif

! Diagnostics follow...
Expand Down
34 changes: 20 additions & 14 deletions src/parameterizations/vertical/MOM_diabatic_aux.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -901,22 +901,24 @@ subroutine applyBoundaryFluxesInOut(CS, G, GV, US, dt, fluxes, optics, nsw, h, t
SurfPressure, & ! Surface pressure (approximated as 0.0) [Pa]
dRhodT, & ! change in density per change in temperature [R degC-1 ~> kg m-3 degC-1]
dRhodS, & ! change in density per change in salinity [R ppt-1 ~> kg m-3 ppt-1]
netheat_rate, & ! netheat but for dt=1 [degC H s-1 ~> degC m s-1 or degC kg m-2 s-1]
netheat_rate, & ! netheat but for dt=1 [degC 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)
! [ppt H s-1 ~> ppt m s-1 or ppt kg m-2 s-1]
netMassInOut_rate! netmassinout but for dt=1 [H s-1 ~> m s-1 or kg m-2 s-1]
! [ppt 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]
real, dimension(SZI_(G), SZK_(G)) :: &
h2d, & ! A 2-d copy of the thicknesses [H ~> m or kg m-2]
T2d, & ! A 2-d copy of the layer temperatures [degC]
pen_TKE_2d, & ! The TKE required to homogenize the heating by shortwave radiation within
! a layer [R Z3 T-2 ~> J m-2]
dSV_dT_2d ! The partial derivative of specific volume with temperature [R-1 degC-1]
real, dimension(SZI_(G),SZK_(G)+1) :: netPen
real, dimension(SZI_(G)) :: &
netPen_rate ! The surface penetrative shortwave heating rate summed over all bands
! [degC H T-1 ~> degC m s-1 or degC kg m-2 s-1]
real, dimension(max(nsw,1),SZI_(G)) :: &
Pen_SW_bnd, & ! The penetrative shortwave heating integrated over a timestep by band
! [degC H ~> degC m or degC kg m-2]
Pen_SW_bnd_rate ! The penetrative shortwave heating rate by band
! [degC H s-1 ~> degC m s-1 or degC kg m-2 s-1]
! [degC H T-1 ~> degC m s-1 or degC kg m-2 s-1]
real, dimension(max(nsw,1),SZI_(G),SZK_(G)) :: &
opacityBand ! The opacity (inverse of the exponential absorption length) of each frequency
! band of shortwave radation in each layer [H-1 ~> m-1 or m2 kg-1]
Expand All @@ -929,7 +931,7 @@ subroutine applyBoundaryFluxesInOut(CS, G, GV, US, dt, fluxes, optics, nsw, h, t
! [Z T-2 R-1 ~> m4 s-2 kg-1]
logical :: calculate_energetics
logical :: calculate_buoyancy
integer :: i, j, is, ie, js, je, k, nz, n
integer :: i, j, is, ie, js, je, k, nz, n, nb
integer :: start, npts
character(len=45) :: mesg

Expand Down Expand Up @@ -970,7 +972,7 @@ 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_in_T, &
!$OMP calculate_buoyancy,netPen,SkinBuoyFlux,GoRho, &
!$OMP calculate_buoyancy,netPen_rate,SkinBuoyFlux,GoRho, &
!$OMP calculate_energetics,dSV_dT,dSV_dS,cTKE,g_Hconv2) &
!$OMP private(opacityBand,h2d,T2d,netMassInOut,netMassOut, &
!$OMP netHeat,netSalt,Pen_SW_bnd,fractionOfForcing, &
Expand Down Expand Up @@ -1334,11 +1336,15 @@ subroutine applyBoundaryFluxesInOut(CS, G, GV, US, dt, fluxes, optics, nsw, h, t
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, US, h2d(:,:), optics_nbands(optics), optics, j, dt_in_T, &
H_limit_fluxes, .true., pen_SW_bnd_rate, netPen)
netPen_rate(:) = 0.0
! Sum over bands and attenuate as a function of depth.
! netPen_rate is the netSW as a function of depth, but only the surface value is used here,
! in which case the values of dt, h, optics and H_limit_fluxes are irrelevant. Consider
! writing a shorter and simpler variant to handle this very limited case.
! call sumSWoverBands(G, GV, US, h2d(:,:), optics_nbands(optics), optics, j, dt_in_T, &
! H_limit_fluxes, .true., pen_SW_bnd_rate, netPen)
do i=is,ie ; do nb=1,nsw ; netPen_rate(i) = netPen_rate(i) + pen_SW_bnd_rate(nb,i) ; enddo ; enddo

! Density derivatives
call calculate_density_derivs(T2d(:,1), tv%S(:,j,1), SurfPressure, &
dRhodT, dRhodS, start, npts, tv%eqn_of_state, scale=US%kg_m3_to_R)
Expand All @@ -1348,9 +1354,9 @@ subroutine applyBoundaryFluxesInOut(CS, G, GV, US, dt, fluxes, optics, nsw, h, t
! 3. Convert to a buoyancy flux, excluding penetrating SW heating
! BGR-Jul 5, 2017: The contribution of SW heating here needs investigated for ePBL.
do i=is,ie
SkinBuoyFlux(i,j) = - GoRho * GV%H_to_Z * US%T_to_s * &
SkinBuoyFlux(i,j) = - GoRho * GV%H_to_Z * &
(dRhodS(i) * (netSalt_rate(i) - tv%S(i,j,1)*netMassInOut_rate(i)) + &
dRhodT(i) * ( netHeat_rate(i) + netPen(i,1)) ) ! [Z2 T-3 ~> m2 s-3]
dRhodT(i) * ( netHeat_rate(i) + netPen_rate(i)) ) ! [Z2 T-3 ~> m2 s-3]
enddo
endif

Expand Down

0 comments on commit 965c2c2

Please sign in to comment.