From 965c2c290a09191258882fe9d4478892ce554c91 Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Sat, 5 Oct 2019 09:38:05 -0400 Subject: [PATCH] +Rescaled _rate variables from extractFluxes1d 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. --- src/core/MOM_forcing_type.F90 | 35 ++++++++++--------- .../vertical/MOM_diabatic_aux.F90 | 34 ++++++++++-------- 2 files changed, 39 insertions(+), 30 deletions(-) diff --git a/src/core/MOM_forcing_type.F90 b/src/core/MOM_forcing_type.F90 index cc94e446cd..f559631606 100644 --- a/src/core/MOM_forcing_type.F90 +++ b/src/core/MOM_forcing_type.F90 @@ -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 @@ -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) @@ -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 @@ -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) ) & @@ -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. @@ -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 @@ -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 @@ -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 @@ -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... diff --git a/src/parameterizations/vertical/MOM_diabatic_aux.F90 b/src/parameterizations/vertical/MOM_diabatic_aux.F90 index ad2f57f2d4..e614524baa 100644 --- a/src/parameterizations/vertical/MOM_diabatic_aux.F90 +++ b/src/parameterizations/vertical/MOM_diabatic_aux.F90 @@ -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] @@ -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 @@ -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, & @@ -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) @@ -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