Skip to content

Commit

Permalink
+Rescale fluxes passed to KPP_NonLocalTransport
Browse files Browse the repository at this point in the history 
  Applied appropriate rescaling for dimensional consistency testing to the net
heat and salt fluxes calculated in calculateBuoyancyFlux and used in the two
KPP_NonLocalTransport_...() routines.  Several comments describing variables
were also either corrected or added.  All answers are bitwise identical.
  • Loading branch information
@Hallberg-NOAA @marshallward
Hallberg-NOAA authored and marshallward committed Nov 28, 2021
1 parent c54dcc6 commit a65fbed
Show file tree
Hide file tree
Showing 4 changed files with 75 additions and 69 deletions.
52 changes: 26 additions & 26 deletions src/core/MOM_forcing_type.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -925,34 +925,32 @@ subroutine calculateBuoyancyFlux1d(G, GV, US, fluxes, optics, nsw, h, Temp, Salt
real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: Salt !< salinity [ppt]
type(thermo_var_ptrs), intent(inout) :: tv !< thermodynamics type
integer, intent(in) :: j !< j-row to work on
real, dimension(SZI_(G),SZK_(GV)+1), intent(inout) :: buoyancyFlux !< buoyancy fluxes [L2 T-3 ~> m2 s-3]
real, dimension(SZI_(G)), intent(inout) :: netHeatMinusSW !< Surface heat flux excluding shortwave
!! [degC H s-1 ~> degC m s-1 or degC kg m-2 s-1]
real, dimension(SZI_(G)), intent(inout) :: netSalt !< surface salt flux
!! [ppt H s-1 ~> ppt m s-1 or ppt kg m-2 s-1]
real, dimension(SZI_(G),SZK_(GV)+1), intent(out) :: buoyancyFlux !< buoyancy fluxes [L2 T-3 ~> m2 s-3]
real, dimension(SZI_(G)), intent(out) :: netHeatMinusSW !< Surface heat flux excluding shortwave
!! [degC H T-1 ~> degC m s-1 or degC kg m-2 s-1]
real, dimension(SZI_(G)), intent(out) :: netSalt !< surface salt flux
!! [ppt H T-1 ~> ppt m s-1 or ppt kg m-2 s-1]

! local variables
integer :: k
real, parameter :: dt = 1. ! to return a rate from extractFluxes1d
real, dimension(SZI_(G)) :: netH ! net FW flux [H s-1 ~> m s-1 or kg m-2 s-1]
real, dimension(SZI_(G)) :: netH ! net FW flux [H T-1 ~> m s-1 or kg m-2 s-1]
real, dimension(SZI_(G)) :: netEvap ! net FW flux leaving ocean via evaporation
! [H s-1 ~> m s-1 or kg m-2 s-1]
real, dimension(SZI_(G)) :: netHeat ! net temp flux [degC H s-1 ~> degC m s-1 or degC kg m-2 s-1]
! [H T-1 ~> m s-1 or kg m-2 s-1]
real, dimension(SZI_(G)) :: netHeat ! net temp flux [degC H T-1 ~> degC m s-1 or degC kg m-2 s-1]
real, dimension(max(nsw,1), SZI_(G)) :: penSWbnd ! penetrating SW radiation 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(SZI_(G)) :: pressure ! pressure at the surface [R L2 T-2 ~> Pa]
real, dimension(SZI_(G)) :: dRhodT ! density partial derivative wrt temp [R degC-1 ~> kg m-3 degC-1]
real, dimension(SZI_(G)) :: dRhodS ! density partial derivative wrt saln [R ppt-1 ~> kg m-3 ppt-1]
real, dimension(SZI_(G),SZK_(GV)+1) :: netPen ! The net penetrating shortwave radiation at each level
! [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 :: useRiverHeatContent
logical :: useCalvingHeatContent
real :: depthBeforeScalingFluxes ! A depth scale [H ~> m or kg m-2]
real :: GoRho ! The gravitational acceleration divided by mean density times some
! unit conversion factors [L2 H-1 s R-1 T-3 ~> m4 kg-1 s-2 or m7 kg-2 s-2]
real :: GoRho ! The gravitational acceleration divided by mean density times a
! unit conversion factor [L2 H-1 R-1 T-2 ~> m4 kg-1 s-2 or m7 kg-2 s-2]
real :: H_limit_fluxes ! Another depth scale [H ~> m or kg m-2]
integer :: i
integer :: i, k

! smg: what do we do when have heat fluxes from calving and river?
useRiverHeatContent = .False.
Expand All @@ -961,38 +959,40 @@ subroutine calculateBuoyancyFlux1d(G, GV, US, fluxes, optics, nsw, h, Temp, Salt
depthBeforeScalingFluxes = max( GV%Angstrom_H, 1.e-30*GV%m_to_H )
pressure(:) = 0.
if (associated(tv%p_surf)) then ; do i=G%isc,G%iec ; pressure(i) = tv%p_surf(i,j) ; enddo ; endif
GoRho = (GV%g_Earth * GV%H_to_Z*US%T_to_s) / GV%Rho0
GoRho = (GV%g_Earth * GV%H_to_Z) / GV%Rho0

H_limit_fluxes = depthBeforeScalingFluxes

! The surface forcing is contained in the fluxes type.
! We aggregate the thermodynamic forcing for a time step into the following:
! netH = water added/removed via surface fluxes [H s-1 ~> m s-1 or kg m-2 s-1]
! netHeat = heat via surface fluxes [degC H s-1 ~> degC m s-1 or degC kg m-2 s-1]
! netSalt = salt via surface fluxes [ppt H s-1 ~> ppt m s-1 or gSalt m-2 s-1]
! netH = water added/removed via surface fluxes [H T-1 ~> m s-1 or kg m-2 s-1]
! netHeat = heat via surface fluxes [degC H T-1 ~> degC m s-1 or degC kg m-2 s-1]
! netSalt = salt via surface fluxes [ppt H T-1 ~> ppt m s-1 or gSalt m-2 s-1]
! Note that unlike other calls to extractFLuxes1d() that return the time-integrated flux
! this call returns the rate because dt=1
call extractFluxes1d(G, GV, US, fluxes, optics, nsw, j, dt*US%s_to_T, &
! this call returns the rate because dt=1 (in arbitrary time units)
call extractFluxes1d(G, GV, US, fluxes, optics, nsw, j, 1.0, &
depthBeforeScalingFluxes, useRiverHeatContent, useCalvingHeatContent, &
h(:,j,:), Temp(:,j,:), netH, netEvap, netHeatMinusSW, &
netSalt, penSWbnd, tv, .false.)

! Sum over bands and attenuate as a function of depth
! netPen is the netSW as a function of depth
call sumSWoverBands(G, GV, US, h(:,j,:), optics_nbands(optics), optics, j, dt*US%s_to_T, &
call sumSWoverBands(G, GV, US, h(:,j,:), optics_nbands(optics), optics, j, 1.0, &
H_limit_fluxes, .true., penSWbnd, netPen)

! Density derivatives
call calculate_density_derivs(Temp(:,j,1), Salt(:,j,1), pressure, dRhodT, dRhodS, &
tv%eqn_of_state, EOS_domain(G%HI))

! Adjust netSalt to reflect dilution effect of FW flux
netSalt(G%isc:G%iec) = netSalt(G%isc:G%iec) - Salt(G%isc:G%iec,j,1) * netH(G%isc:G%iec) ! ppt H/s
! [ppt H T-1 ~> ppt m s-1 or ppt kg m-2 s-1]
netSalt(G%isc:G%iec) = netSalt(G%isc:G%iec) - Salt(G%isc:G%iec,j,1) * netH(G%isc:G%iec)

! Add in the SW heating for purposes of calculating the net
! surface buoyancy flux affecting the top layer.
! [degC H T-1 ~> degC m s-1 or degC kg m-2 s-1]
!netHeat(:) = netHeatMinusSW(:) + sum( penSWbnd, dim=1 )
netHeat(G%isc:G%iec) = netHeatMinusSW(G%isc:G%iec) + netPen(G%isc:G%iec,1) ! K H/s
netHeat(G%isc:G%iec) = netHeatMinusSW(G%isc:G%iec) + netPen(G%isc:G%iec,1)

! Convert to a buoyancy flux, excluding penetrating SW heating
buoyancyFlux(G%isc:G%iec,1) = - GoRho * ( dRhodS(G%isc:G%iec) * netSalt(G%isc:G%iec) + &
Expand Down Expand Up @@ -1020,9 +1020,9 @@ subroutine calculateBuoyancyFlux2d(G, GV, US, fluxes, optics, h, Temp, Salt, tv,
type(thermo_var_ptrs), intent(inout) :: tv !< thermodynamics type
real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1), intent(inout) :: buoyancyFlux !< buoyancy fluxes [L2 T-3 ~> m2 s-3]
real, dimension(SZI_(G),SZJ_(G)), intent(inout) :: netHeatMinusSW !< surface heat flux excluding shortwave
!! [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),SZJ_(G)), intent(inout) :: netSalt !< Net surface salt flux
!! [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]
! local variables
integer :: j

Expand Down
57 changes: 32 additions & 25 deletions src/parameterizations/vertical/MOM_CVMix_KPP.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -538,9 +538,9 @@ logical function KPP_init(paramFile, G, GV, US, diag, Time, CS, passive)
CS%id_buoyFlux = register_diag_field('ocean_model', 'KPP_buoyFlux', diag%axesTi, Time, &
'Surface (and penetrating) buoyancy flux, as used by [CVMix] KPP', 'm2/s3', conversion=US%L_to_m**2*US%s_to_T**3)
CS%id_QminusSW = register_diag_field('ocean_model', 'KPP_QminusSW', diag%axesT1, Time, &
'Net temperature flux ignoring short-wave, as used by [CVMix] KPP', 'K m/s')
'Net temperature flux ignoring short-wave, as used by [CVMix] KPP', 'K m/s', conversion=GV%H_to_m*US%s_to_T)
CS%id_netS = register_diag_field('ocean_model', 'KPP_netSalt', diag%axesT1, Time, &
'Effective net surface salt flux, as used by [CVMix] KPP', 'ppt m/s')
'Effective net surface salt flux, as used by [CVMix] KPP', 'ppt m/s', conversion=GV%H_to_m*US%s_to_T)
CS%id_Kt_KPP = register_diag_field('ocean_model', 'KPP_Kheat', diag%axesTi, Time, &
'Heat diffusivity due to KPP, as calculated by [CVMix] KPP', 'm2/s')
CS%id_Kd_in = register_diag_field('ocean_model', 'KPP_Kd_in', diag%axesTi, Time, &
Expand All @@ -554,13 +554,17 @@ logical function KPP_init(paramFile, G, GV, US, diag, Time, CS, passive)
CS%id_NLTs = register_diag_field('ocean_model', 'KPP_NLtransport_salt', diag%axesTi, Time, &
'Non-local tranpsort (Cs*G(sigma)) for scalars, as calculated by [CVMix] KPP', 'nondim')
CS%id_NLT_dTdt = register_diag_field('ocean_model', 'KPP_NLT_dTdt', diag%axesTL, Time, &
'Temperature tendency due to non-local transport of heat, as calculated by [CVMix] KPP', 'K/s')
'Temperature tendency due to non-local transport of heat, as calculated by [CVMix] KPP', &
'K/s', conversion=US%s_to_T)
CS%id_NLT_dSdt = register_diag_field('ocean_model', 'KPP_NLT_dSdt', diag%axesTL, Time, &
'Salinity tendency due to non-local transport of salt, as calculated by [CVMix] KPP', 'ppt/s')
'Salinity tendency due to non-local transport of salt, as calculated by [CVMix] KPP', &
'ppt/s', conversion=US%s_to_T)
CS%id_NLT_temp_budget = register_diag_field('ocean_model', 'KPP_NLT_temp_budget', diag%axesTL, Time, &
'Heat content change due to non-local transport, as calculated by [CVMix] KPP', 'W/m^2')
'Heat content change due to non-local transport, as calculated by [CVMix] KPP', &
'W/m^2', conversion=US%QRZ_T_to_W_m2)
CS%id_NLT_saln_budget = register_diag_field('ocean_model', 'KPP_NLT_saln_budget', diag%axesTL, Time, &
'Salt content change due to non-local transport, as calculated by [CVMix] KPP', 'kg/(sec*m^2)')
'Salt content change due to non-local transport, as calculated by [CVMix] KPP', &
'kg/(sec*m^2)', conversion=US%RZ_T_to_kg_m2s)
CS%id_Tsurf = register_diag_field('ocean_model', 'KPP_Tsurf', diag%axesT1, Time, &
'Temperature of surface layer (10% of OBL depth) as passed to [CVMix] KPP', 'C')
CS%id_Ssurf = register_diag_field('ocean_model', 'KPP_Ssurf', diag%axesT1, Time, &
Expand Down Expand Up @@ -1179,7 +1183,7 @@ subroutine KPP_compute_BLD(CS, G, GV, US, h, Temp, Salt, u, v, tv, uStar, buoyFl
! Calculate Bulk Richardson number from eq (21) of LMD94
BulkRi_1d = CVmix_kpp_compute_bulk_Richardson( &
zt_cntr = cellHeight(1:GV%ke), & ! Depth of cell center [m]
delta_buoy_cntr=GoRho*deltaRho, & ! Bulk buoyancy difference, Br-B(z) [s-1]
delta_buoy_cntr=GoRho*deltaRho, & ! Bulk buoyancy difference, Br-B(z) [m s-2]
delta_Vsqr_cntr=deltaU2, & ! Square of resolved velocity difference [m2 s-2]
ws_cntr=Ws_1d, & ! Turbulent velocity scale profile [m s-1]
N_iface=CS%N(i,j,:), & ! Buoyancy frequency [s-1]
Expand Down Expand Up @@ -1285,12 +1289,12 @@ subroutine KPP_smooth_BLD(CS,G,GV,h)
real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h !< Layer/level thicknesses [H ~> m or kg m-2]

! local
real, dimension(SZI_(G),SZJ_(G)) :: OBLdepth_prev ! OBLdepth before s.th smoothing iteration
real, dimension(SZI_(G),SZJ_(G)) :: OBLdepth_prev ! OBLdepth before s.th smoothing iteration [m]
real, dimension( GV%ke ) :: cellHeight ! Cell center heights referenced to surface [m]
! (negative in the ocean)
real, dimension( GV%ke+1 ) :: iFaceHeight ! Interface heights referenced to surface [m]
! (negative in the ocean)
real :: wc, ww, we, wn, ws ! averaging weights for smoothing
real :: wc, ww, we, wn, ws ! averaging weights for smoothing [nondim]
real :: dh ! The local thickness used for calculating interface positions [m]
real :: hcorr ! A cumulative correction arising from inflation of vanished layers [m]
integer :: i, j, k, s
Expand Down Expand Up @@ -1390,13 +1394,14 @@ subroutine KPP_NonLocalTransport_temp(CS, G, GV, h, nonLocalTrans, surfFlux, &
real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h !< Layer/level thickness [H ~> m or kg m-2]
real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1), intent(in) :: nonLocalTrans !< Non-local transport [nondim]
real, dimension(SZI_(G),SZJ_(G)), intent(in) :: surfFlux !< Surface flux of temperature
!! [degC H s-1 ~> degC m s-1 or degC kg m-2 s-1]
real, intent(in) :: dt !< Time-step [s]
real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(inout) :: scalar !< temperature
real, intent(in) :: C_p !< Seawater specific heat capacity [J kg-1 degC-1]
!! [degC H T-1 ~> degC m s-1 or degC kg m-2 s-1]
real, intent(in) :: dt !< Time-step [T ~> s]
real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(inout) :: scalar !< temperature [degC]
real, intent(in) :: C_p !< Seawater specific heat capacity
!! [Q degC-1 ~> J kg-1 degC-1]

integer :: i, j, k
real, dimension( SZI_(G), SZJ_(G),SZK_(GV) ) :: dtracer
real, dimension( SZI_(G), SZJ_(G),SZK_(GV) ) :: dtracer ! Rate of tracer change [degC T-1 ~> degC s-1]


dtracer(:,:,:) = 0.0
Expand Down Expand Up @@ -1431,8 +1436,9 @@ subroutine KPP_NonLocalTransport_temp(CS, G, GV, h, nonLocalTrans, surfFlux, &
do k = 1, GV%ke
do j = G%jsc, G%jec
do i = G%isc, G%iec
! Here dtracer has units of [Q R Z T-1 ~> W m-2].
dtracer(i,j,k) = (nonLocalTrans(i,j,k) - nonLocalTrans(i,j,k+1)) * &
surfFlux(i,j) * C_p * GV%H_to_kg_m2
surfFlux(i,j) * C_p * GV%H_to_RZ
enddo
enddo
enddo
Expand All @@ -1446,18 +1452,18 @@ end subroutine KPP_NonLocalTransport_temp
!> This routine is a useful prototype for other material tracers.
subroutine KPP_NonLocalTransport_saln(CS, G, GV, h, nonLocalTrans, surfFlux, dt, scalar)

type(KPP_CS), intent(in) :: CS !< Control structure
type(ocean_grid_type), intent(in) :: G !< Ocean grid
type(verticalGrid_type), intent(in) :: GV !< Ocean vertical grid
real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h !< Layer/level thickness [H ~> m or kg m-2]
type(KPP_CS), intent(in) :: CS !< Control structure
type(ocean_grid_type), intent(in) :: G !< Ocean grid
type(verticalGrid_type), intent(in) :: GV !< Ocean vertical grid
real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h !< Layer/level thickness [H ~> m or kg m-2]
real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1), intent(in) :: nonLocalTrans !< Non-local transport [nondim]
real, dimension(SZI_(G),SZJ_(G)), intent(in) :: surfFlux !< Surface flux of salt
!! [ppt H s-1 ~> ppt m s-1 or ppt kg m-2 s-1]
real, intent(in) :: dt !< Time-step [s]
real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(inout) :: scalar !< Scalar (scalar units [conc])
real, dimension(SZI_(G),SZJ_(G)), intent(in) :: surfFlux !< Surface flux of salt
!! [ppt H T-1 ~> ppt m s-1 or ppt kg m-2 s-1]
real, intent(in) :: dt !< Time-step [T ~> s]
real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(inout) :: scalar !< Salinity [ppt]

integer :: i, j, k
real, dimension( SZI_(G), SZJ_(G),SZK_(GV) ) :: dtracer
real, dimension( SZI_(G), SZJ_(G),SZK_(GV) ) :: dtracer ! Rate of tracer change [ppt T-1 ~> ppt s-1]


dtracer(:,:,:) = 0.0
Expand Down Expand Up @@ -1492,8 +1498,9 @@ subroutine KPP_NonLocalTransport_saln(CS, G, GV, h, nonLocalTrans, surfFlux, dt,
do k = 1, GV%ke
do j = G%jsc, G%jec
do i = G%isc, G%iec
! Here dtracer has units of [ppt R Z T-1 ~> ppt kg m-2 s-1]
dtracer(i,j,k) = (nonLocalTrans(i,j,k) - nonLocalTrans(i,j,k+1)) * &
surfFlux(i,j) * GV%H_to_kg_m2
surfFlux(i,j) * GV%H_to_RZ
enddo
enddo
enddo
Expand Down
Loading

0 comments on commit a65fbed

Please sign in to comment.