Skip to content

Commit

Permalink
Rescaled density units in MOM_mixed_layer_restrat
Browse files Browse the repository at this point in the history 
  Rescaled density units in MOM_mixed_layer_restrat for dimensional consistency
testing.  All answers are bitwise identical.
  • Loading branch information
@Hallberg-NOAA
Hallberg-NOAA committed Sep 30, 2019
1 parent 3fb51fd commit 04901cd
Showing 1 changed file with 18 additions and 15 deletions.
33 changes: 18 additions & 15 deletions src/parameterizations/lateral/MOM_mixed_layer_restrat.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -50,7 +50,7 @@ module MOM_mixed_layer_restrat
!! based on the parameter MLE_DENSITY_DIFF.
real :: MLE_MLD_decay_time !< Time-scale to use in a running-mean when MLD is retreating [T ~> s].
real :: MLE_MLD_decay_time2 !< Time-scale to use in a running-mean when filtered MLD is retreating [T ~> s].
real :: MLE_density_diff !< Density difference used in detecting mixed-layer depth [kg m-3].
real :: MLE_density_diff !< Density difference used in detecting mixed-layer depth [R ~> kg m-3].
real :: MLE_tail_dh !< Fraction by which to extend the mixed-layer restratification
!! depth used for a smoother stream function at the base of
!! the mixed-layer [nondim].
Expand Down Expand Up @@ -147,8 +147,8 @@ subroutine mixedlayer_restrat_general(h, uhtr, vhtr, tv, forces, dt_in_T, MLD_in
MLD_slow, & ! Mixed layer depth actually used in MLE restratification parameterization [H ~> m or kg m-2]
htot_slow, & ! The sum of the thicknesses of layers in the mixed layer [H ~> m or kg m-2]
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 [m3 L2 Z-1 T-2 kg-1 ~> m4 s-2 kg-1]
real :: rho_ml(SZI_(G)) ! Potential density relative to the surface [kg m-3]
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 :: h_vel ! htot interpolated onto velocity points [Z ~> m] (not H).
Expand All @@ -174,11 +174,12 @@ subroutine mixedlayer_restrat_general(h, uhtr, vhtr, tv, forces, dt_in_T, MLD_in
! for diagnostic purposes.
real :: uDml_diag(SZIB_(G),SZJ_(G)), vDml_diag(SZI_(G),SZJB_(G))
integer :: i, j, k, is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz
real, dimension(SZI_(G)) :: rhoSurf, deltaRhoAtKm1, deltaRhoAtK
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)) :: rhoAtK, rho1, d1, pRef_N2 ! Used for N2
real :: aFac, bFac, ddRho
real :: aFac, bFac ! Nondimensional ratios [nondim]
real :: ddRho ! A density difference [R ~> kg m-3]
real :: hAtVel, zpa, zpb, dh, res_scaling_fac
real :: I_LFront ! The inverse of the frontal length scale [L-1 ~> m-1]
logical :: proper_averaging, line_is_empty, keep_going, res_upscale
Expand All @@ -205,15 +206,17 @@ subroutine mixedlayer_restrat_general(h, uhtr, vhtr, tv, forces, dt_in_T, MLD_in
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)
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)
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)
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)
deltaRhoAtK(:) = deltaRhoAtK(:) - rhoSurf(:) ! Density difference between layer K and surface
do i = is-1, ie+1
ddRho = deltaRhoAtK(i) - deltaRhoAtKm1(i)
Expand Down Expand Up @@ -282,7 +285,7 @@ subroutine mixedlayer_restrat_general(h, uhtr, vhtr, tv, forces, dt_in_T, MLD_in
uDml(:) = 0.0 ; vDml(:) = 0.0
uDml_slow(:) = 0.0 ; vDml_slow(:) = 0.0
I4dt = 0.25 / (dt_in_T)
g_Rho0 = GV%g_Earth / (US%R_to_kg_m3*GV%Rho0)
g_Rho0 = GV%g_Earth / GV%Rho0
h_neglect = GV%H_subroundoff
dz_neglect = GV%H_subroundoff*GV%H_to_Z
proper_averaging = .not. CS%MLE_use_MLD_ave_bug
Expand Down Expand Up @@ -316,7 +319,7 @@ subroutine mixedlayer_restrat_general(h, uhtr, vhtr, tv, forces, dt_in_T, MLD_in
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)
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)
line_is_empty = .true.
do i=is-1,ie+1
if (htot_fast(i,j) < MLD_fast(i,j)) then
Expand Down Expand Up @@ -578,8 +581,8 @@ subroutine mixedlayer_restrat_BML(h, uhtr, vhtr, tv, forces, dt_in_T, G, GV, US,
real, dimension(SZI_(G),SZJ_(G)) :: &
htot, & ! The sum of the thicknesses of layers in the mixed layer [H ~> m or kg m-2]
Rml_av ! g_Rho0 times the average mixed layer density [L2 Z-1 T-2 ~> m s-2]
real :: g_Rho0 ! G_Earth/Rho0 [m3 L2 Z-1 T-2 kg-1 ~> m4 s-2 kg-1]
real :: Rho0(SZI_(G)) ! Potential density relative to the surface [kg m-3]
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 :: h_vel ! htot interpolated onto velocity points [Z ~> m]. (The units are not H.)
Expand Down Expand Up @@ -616,7 +619,7 @@ subroutine mixedlayer_restrat_BML(h, uhtr, vhtr, tv, forces, dt_in_T, G, GV, US,

uDml(:) = 0.0 ; vDml(:) = 0.0
I4dt = 0.25 / (dt_in_T)
g_Rho0 = GV%g_Earth / (US%R_to_kg_m3*GV%Rho0)
g_Rho0 = GV%g_Earth / GV%Rho0
use_EOS = associated(tv%eqn_of_state)
h_neglect = GV%H_subroundoff
dz_neglect = GV%H_subroundoff*GV%H_to_Z
Expand All @@ -640,7 +643,7 @@ subroutine mixedlayer_restrat_BML(h, uhtr, vhtr, tv, forces, dt_in_T, G, GV, US,
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)
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)
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 Expand Up @@ -821,7 +824,7 @@ logical function mixedlayer_restrat_init(Time, G, GV, US, param_file, diag, CS,

! Nonsense values to cause problems when these parameters are not used
CS%MLE_MLD_decay_time = -9.e9*US%s_to_T
CS%MLE_density_diff = -9.e9
CS%MLE_density_diff = -9.e9*US%kg_m3_to_R
CS%MLE_tail_dh = -9.e9
CS%MLE_use_PBL_MLD = .false.
CS%MLE_MLD_stretch = -9.e9
Expand Down Expand Up @@ -867,7 +870,7 @@ logical function mixedlayer_restrat_init(Time, G, GV, US, param_file, diag, CS,
call get_param(param_file, mdl, "MLE_DENSITY_DIFF", CS%MLE_density_diff, &
"Density difference used to detect the mixed-layer "//&
"depth used for the mixed-layer eddy parameterization "//&
"by Fox-Kemper et al. (2010)", units="kg/m3", default=0.03)
"by Fox-Kemper et al. (2010)", units="kg/m3", default=0.03, scale=US%kg_m3_to_R)
endif
call get_param(param_file, mdl, "MLE_TAIL_DH", CS%MLE_tail_dh, &
"Fraction by which to extend the mixed-layer restratification "//&
Expand Down

0 comments on commit 04901cd

Please sign in to comment.