Skip to content

Commit

Permalink
Rescaled pressures in 17 calculate_density calls
Browse files Browse the repository at this point in the history 
  Rescaled pressures in 17 calculate_density or related calls and internal
pressure variables in MOM_diabatic_aux.  All answers are bitwise identical.
  • Loading branch information
@Hallberg-NOAA
Hallberg-NOAA committed Apr 13, 2020
1 parent 8df34f4 commit 500c9e3
Show file tree
Hide file tree
Showing 6 changed files with 62 additions and 73 deletions.
19 changes: 7 additions & 12 deletions src/core/MOM_forcing_type.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -851,10 +851,7 @@ subroutine extractFluxes2d(G, GV, US, fluxes, optics, nsw, dt, FluxRescaleDepth,
logical, intent(in) :: aggregate_FW !< For determining how to aggregate the forcing.

integer :: j
!$OMP parallel do default(none) shared(G, GV, US, fluxes, optics, nsw, dt, FluxRescaleDepth, &
!$OMP useRiverHeatContent, useCalvingHeatContent, &
!$OMP h,T,netMassInOut,netMassOut,Net_heat,Net_salt,Pen_SW_bnd,tv, &
!$OMP aggregate_FW)
!$OMP parallel do default(shared)
do j=G%jsc, G%jec
call extractFluxes1d(G, GV, US, fluxes, optics, nsw, j, dt, &
FluxRescaleDepth, useRiverHeatContent, useCalvingHeatContent,&
Expand Down Expand Up @@ -891,15 +888,15 @@ subroutine calculateBuoyancyFlux1d(G, GV, US, fluxes, optics, nsw, h, Temp, Salt
logical, optional, intent(in) :: skip_diags !< If present and true, skip calculating
!! diagnostics inside extractFluxes1d()
! local variables
integer :: start, npts, k
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)) :: 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-2 or degC kg m-2 s-1]
real, dimension(max(nsw,1), SZI_(G)) :: penSWbnd ! penetrating SW radiation by band
! [degC H ~> degC m or degC kg m-2]
real, dimension(SZI_(G)) :: pressure ! pressurea the surface [Pa]
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_(G)+1) :: netPen ! The net penetrating shortwave radiation at each level
Expand All @@ -917,10 +914,8 @@ subroutine calculateBuoyancyFlux1d(G, GV, US, fluxes, optics, nsw, h, Temp, Salt
useCalvingHeatContent = .False.

depthBeforeScalingFluxes = max( GV%Angstrom_H, 1.e-30*GV%m_to_H )
pressure(:) = 0. ! Ignore atmospheric pressure
pressure(:) = 0. ! Ignores atmospheric pressure ###
GoRho = (GV%g_Earth * GV%H_to_Z*US%T_to_s) / GV%Rho0
start = 1 + G%isc - G%isd
npts = 1 + G%iec - G%isc

H_limit_fluxes = depthBeforeScalingFluxes

Expand All @@ -931,7 +926,7 @@ subroutine calculateBuoyancyFlux1d(G, GV, US, fluxes, optics, nsw, h, Temp, Salt
! netSalt = salt via surface fluxes [ppt H s-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, &
call extractFluxes1d(G, GV, US, fluxes, optics, nsw, j, dt*US%s_to_T, &
depthBeforeScalingFluxes, useRiverHeatContent, useCalvingHeatContent, &
h(:,j,:), Temp(:,j,:), netH, netEvap, netHeatMinusSW, &
netSalt, penSWbnd, tv, .false., skip_diags=skip_diags)
Expand All @@ -942,8 +937,8 @@ subroutine calculateBuoyancyFlux1d(G, GV, US, fluxes, optics, nsw, h, Temp, Salt
H_limit_fluxes, .true., penSWbnd, netPen)

! Density derivatives
call calculate_density_derivs(Temp(:,j,1), Salt(:,j,1), pressure, &
dRhodT, dRhodS, start, npts, tv%eqn_of_state, scale=US%kg_m3_to_R)
call calculate_density_derivs(Temp(:,j,1), Salt(:,j,1), pressure, dRhodT, dRhodS, G%HI, &
tv%eqn_of_state, US)

! 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
Expand Down
19 changes: 10 additions & 9 deletions src/parameterizations/lateral/MOM_mixed_layer_restrat.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -149,7 +149,7 @@ subroutine mixedlayer_restrat_general(h, uhtr, vhtr, tv, forces, dt, MLD_in, Var
Rml_av_slow ! g_Rho0 times the average mixed layer density [L2 Z-1 T-2 ~> m s-2]
real :: g_Rho0 ! G_Earth/Rho0 [L2 Z-1 T-2 R-1 ~> m4 s-2 kg-1]
real :: rho_ml(SZI_(G)) ! Potential density relative to the surface [R ~> kg m-3]
real :: p0(SZI_(G)) ! A pressure of 0 [Pa]
real :: p0(SZI_(G)) ! A pressure of 0 [R L2 T-2 ~> Pa]

real :: h_vel ! htot interpolated onto velocity points [Z ~> m] (not H).
real :: absf ! absolute value of f, interpolated to velocity points [T-1 ~> s-1]
Expand All @@ -176,7 +176,8 @@ subroutine mixedlayer_restrat_general(h, uhtr, vhtr, tv, forces, dt, MLD_in, Var
integer :: i, j, k, is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz
real, dimension(SZI_(G)) :: rhoSurf, deltaRhoAtKm1, deltaRhoAtK ! Densities [R ~> kg m-3]
real, dimension(SZI_(G)) :: dK, dKm1 ! Depths of layer centers [H ~> m or kg m-2].
real, dimension(SZI_(G)) :: pRef_MLD ! A reference pressure for calculating the mixed layer densities [Pa].
real, dimension(SZI_(G)) :: pRef_MLD ! A reference pressure for calculating the mixed layer
! densities [R L2 T-2 ~> Pa].
real, dimension(SZI_(G)) :: rhoAtK, rho1, d1, pRef_N2 ! Used for N2
real :: aFac, bFac ! Nondimensional ratios [nondim]
real :: ddRho ! A density difference [R ~> kg m-3]
Expand Down Expand Up @@ -206,17 +207,17 @@ subroutine mixedlayer_restrat_general(h, uhtr, vhtr, tv, forces, dt, MLD_in, Var
pRef_MLD(:) = 0.
do j = js-1, je+1
dK(:) = 0.5 * h(:,j,1) ! Depth of center of surface layer
call calculate_density(tv%T(:,j,1), tv%S(:,j,1), pRef_MLD, rhoSurf, is-1, ie-is+3, &
tv%eqn_of_state, scale=US%kg_m3_to_R)
call calculate_density(tv%T(:,j,1), tv%S(:,j,1), pRef_MLD, rhoSurf, G%HI, &
tv%eqn_of_state, US, halo=1)
deltaRhoAtK(:) = 0.
MLD_fast(:,j) = 0.
do k = 2, nz
dKm1(:) = dK(:) ! Depth of center of layer K-1
dK(:) = dK(:) + 0.5 * ( h(:,j,k) + h(:,j,k-1) ) ! Depth of center of layer K
! Mixed-layer depth, using sigma-0 (surface reference pressure)
deltaRhoAtKm1(:) = deltaRhoAtK(:) ! Store value from previous iteration of K
call calculate_density(tv%T(:,j,k), tv%S(:,j,k), pRef_MLD, deltaRhoAtK, is-1, ie-is+3, &
tv%eqn_of_state, scale=US%kg_m3_to_R)
call calculate_density(tv%T(:,j,k), tv%S(:,j,k), pRef_MLD, deltaRhoAtK, G%HI, &
tv%eqn_of_state, US, halo=1)
do i = is-1,ie+1
deltaRhoAtK(i) = deltaRhoAtK(i) - rhoSurf(i) ! Density difference between layer K and surface
enddo
Expand Down Expand Up @@ -321,7 +322,7 @@ subroutine mixedlayer_restrat_general(h, uhtr, vhtr, tv, forces, dt, MLD_in, Var
h_avail(i,j,k) = max(I4dt*G%areaT(i,j)*(h(i,j,k)-GV%Angstrom_H),0.0)
enddo
if (keep_going) then
call calculate_density(tv%T(:,j,k),tv%S(:,j,k),p0,rho_ml(:),is-1,ie-is+3,tv%eqn_of_state, scale=US%kg_m3_to_R)
call calculate_density(tv%T(:,j,k), tv%S(:,j,k), p0, rho_ml(:), G%HI, tv%eqn_of_state, US, halo=1)
line_is_empty = .true.
do i=is-1,ie+1
if (htot_fast(i,j) < MLD_fast(i,j)) then
Expand Down Expand Up @@ -585,7 +586,7 @@ subroutine mixedlayer_restrat_BML(h, uhtr, vhtr, tv, forces, dt, G, GV, US, CS)
Rml_av ! g_Rho0 times the average mixed layer density [L2 Z-1 T-2 ~> m s-2]
real :: g_Rho0 ! G_Earth/Rho0 [L2 Z-1 T-2 R-1 ~> m4 s-2 kg-1]
real :: Rho0(SZI_(G)) ! Potential density relative to the surface [R ~> kg m-3]
real :: p0(SZI_(G)) ! A pressure of 0 [Pa]
real :: p0(SZI_(G)) ! A pressure of 0 [R L2 T-2 ~> Pa]

real :: h_vel ! htot interpolated onto velocity points [Z ~> m]. (The units are not H.)
real :: absf ! absolute value of f, interpolated to velocity points [T-1 ~> s-1]
Expand Down Expand Up @@ -645,7 +646,7 @@ subroutine mixedlayer_restrat_BML(h, uhtr, vhtr, tv, forces, dt, G, GV, US, CS)
htot(i,j) = 0.0 ; Rml_av(i,j) = 0.0
enddo
do k=1,nkml
call calculate_density(tv%T(:,j,k),tv%S(:,j,k),p0,Rho0(:),is-1,ie-is+3,tv%eqn_of_state, scale=US%kg_m3_to_R)
call calculate_density(tv%T(:,j,k), tv%S(:,j,k), p0, Rho0(:), G%HI, tv%eqn_of_state, US, halo=1)
do i=is-1,ie+1
Rml_av(i,j) = Rml_av(i,j) + h(i,j,k)*Rho0(i)
htot(i,j) = htot(i,j) + h(i,j,k)
Expand Down
24 changes: 12 additions & 12 deletions src/parameterizations/lateral/MOM_thickness_diffuse.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -594,7 +594,7 @@ subroutine thickness_diffuse_full(h, e, Kh_u, Kh_v, tv, uhD, vhD, cg1, dt, G, GV
hN2_x_PE ! thickness in m times Brunt-Vaisala freqeuncy at u-points [L2 Z-1 T-2 ~> m s-2],
! used for calculating PE release
real, dimension(SZI_(G), SZJ_(G), SZK_(G)+1) :: &
pres, & ! The pressure at an interface [Pa].
pres, & ! The pressure at an interface [R L2 T-2 ~> Pa].
h_avail_rsum ! The running sum of h_avail above an interface [H L2 T-1 ~> m3 s-1 or kg s-1].
real, dimension(SZIB_(G)) :: &
drho_dT_u, & ! The derivative of density with temperature at u points [R degC-1 ~> kg m-3 degC-1]
Expand All @@ -607,11 +607,11 @@ subroutine thickness_diffuse_full(h, e, Kh_u, Kh_v, tv, uhD, vhD, cg1, dt, G, GV
real, dimension(SZIB_(G)) :: &
T_u, & ! Temperature on the interface at the u-point [degC].
S_u, & ! Salinity on the interface at the u-point [ppt].
pres_u ! Pressure on the interface at the u-point [Pa].
pres_u ! Pressure on the interface at the u-point [R L2 T-2 ~> Pa].
real, dimension(SZI_(G)) :: &
T_v, & ! Temperature on the interface at the v-point [degC].
S_v, & ! Salinity on the interface at the v-point [ppt].
pres_v ! Pressure on the interface at the v-point [Pa].
pres_v ! Pressure on the interface at the v-point [R L2 T-2 ~> Pa].
real :: Work_u(SZIB_(G), SZJ_(G)) ! The work being done by the thickness
real :: Work_v(SZI_(G), SZJB_(G)) ! diffusion integrated over a cell [R Z L4 T-3 ~> W ]
real :: Work_h ! The work averaged over an h-cell [R Z L2 T-3 ~> W m-2].
Expand Down Expand Up @@ -720,7 +720,7 @@ subroutine thickness_diffuse_full(h, e, Kh_u, Kh_v, tv, uhD, vhD, cg1, dt, G, GV
h_avail(i,j,1) = max(I4dt*G%areaT(i,j)*(h(i,j,1)-GV%Angstrom_H),0.0)
h_avail_rsum(i,j,2) = h_avail(i,j,1)
h_frac(i,j,1) = 1.0
pres(i,j,2) = pres(i,j,1) + GV%H_to_Pa*h(i,j,1)
pres(i,j,2) = pres(i,j,1) + (GV%g_Earth*GV%H_to_RZ) * h(i,j,1)
enddo ; enddo
!$OMP do
do j=js-1,je+1
Expand All @@ -729,7 +729,7 @@ subroutine thickness_diffuse_full(h, e, Kh_u, Kh_v, tv, uhD, vhD, cg1, dt, G, GV
h_avail_rsum(i,j,k+1) = h_avail_rsum(i,j,k) + h_avail(i,j,k)
h_frac(i,j,k) = 0.0 ; if (h_avail(i,j,k) > 0.0) &
h_frac(i,j,k) = h_avail(i,j,k) / h_avail_rsum(i,j,k+1)
pres(i,j,K+1) = pres(i,j,K) + GV%H_to_Pa*h(i,j,k)
pres(i,j,K+1) = pres(i,j,K) + (GV%g_Earth*GV%H_to_RZ) * h(i,j,k)
enddo ; enddo
enddo
!$OMP do
Expand Down Expand Up @@ -778,7 +778,7 @@ subroutine thickness_diffuse_full(h, e, Kh_u, Kh_v, tv, uhD, vhD, cg1, dt, G, GV
S_u(I) = 0.25*((S(i,j,k) + S(i+1,j,k)) + (S(i,j,k-1) + S(i+1,j,k-1)))
enddo
call calculate_density_derivs(T_u, S_u, pres_u, drho_dT_u, &
drho_dS_u, (is-IsdB+1)-1, ie-is+2, tv%eqn_of_state, scale=US%kg_m3_to_R)
drho_dS_u, (is-IsdB+1)-1, ie-is+2, tv%eqn_of_state, US=US)
endif

do I=is-1,ie
Expand Down Expand Up @@ -1028,8 +1028,8 @@ subroutine thickness_diffuse_full(h, e, Kh_u, Kh_v, tv, uhD, vhD, cg1, dt, G, GV
T_v(i) = 0.25*((T(i,j,k) + T(i,j+1,k)) + (T(i,j,k-1) + T(i,j+1,k-1)))
S_v(i) = 0.25*((S(i,j,k) + S(i,j+1,k)) + (S(i,j,k-1) + S(i,j+1,k-1)))
enddo
call calculate_density_derivs(T_v, S_v, pres_v, drho_dT_v, &
drho_dS_v, is, ie-is+1, tv%eqn_of_state, scale=US%kg_m3_to_R)
call calculate_density_derivs(T_v, S_v, pres_v, drho_dT_v, drho_dS_v, G%HI, &
tv%eqn_of_state, US)
endif
do i=is,ie
if (calc_derivatives) then
Expand Down Expand Up @@ -1260,8 +1260,8 @@ subroutine thickness_diffuse_full(h, e, Kh_u, Kh_v, tv, uhD, vhD, cg1, dt, G, GV
T_u(I) = 0.5*(T(i,j,1) + T(i+1,j,1))
S_u(I) = 0.5*(S(i,j,1) + S(i+1,j,1))
enddo
call calculate_density_derivs(T_u, S_u, pres_u, drho_dT_u, &
drho_dS_u, (is-IsdB+1)-1, ie-is+2, tv%eqn_of_state, scale=US%kg_m3_to_R)
call calculate_density_derivs(T_u, S_u, pres_u, drho_dT_u, drho_dS_u, &
(is-IsdB+1)-1, ie-is+2, tv%eqn_of_state, US=US)
endif
do I=is-1,ie
uhD(I,j,1) = -uhtot(I,j)
Expand Down Expand Up @@ -1290,8 +1290,8 @@ subroutine thickness_diffuse_full(h, e, Kh_u, Kh_v, tv, uhD, vhD, cg1, dt, G, GV
T_v(i) = 0.5*(T(i,j,1) + T(i,j+1,1))
S_v(i) = 0.5*(S(i,j,1) + S(i,j+1,1))
enddo
call calculate_density_derivs(T_v, S_v, pres_v, drho_dT_v, &
drho_dS_v, is, ie-is+1, tv%eqn_of_state, scale=US%kg_m3_to_R)
call calculate_density_derivs(T_v, S_v, pres_v, drho_dT_v, drho_dS_v, G%HI, &
tv%eqn_of_state, US)
endif
do i=is,ie
vhD(i,J,1) = -vhtot(i,J)
Expand Down
Loading

0 comments on commit 500c9e3

Please sign in to comment.