From 04901cd0ae3a14c5d1e8397b5be507d8d2413dd8 Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Mon, 30 Sep 2019 07:25:36 -0400 Subject: [PATCH] Rescaled density units in MOM_mixed_layer_restrat Rescaled density units in MOM_mixed_layer_restrat for dimensional consistency testing. All answers are bitwise identical. --- .../lateral/MOM_mixed_layer_restrat.F90 | 33 ++++++++++--------- 1 file changed, 18 insertions(+), 15 deletions(-) diff --git a/src/parameterizations/lateral/MOM_mixed_layer_restrat.F90 b/src/parameterizations/lateral/MOM_mixed_layer_restrat.F90 index ca62160bc1..182ec2dc0c 100644 --- a/src/parameterizations/lateral/MOM_mixed_layer_restrat.F90 +++ b/src/parameterizations/lateral/MOM_mixed_layer_restrat.F90 @@ -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]. @@ -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). @@ -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 @@ -205,7 +206,8 @@ 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 @@ -213,7 +215,8 @@ subroutine mixedlayer_restrat_general(h, uhtr, vhtr, tv, forces, dt_in_T, MLD_in 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) @@ -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 @@ -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 @@ -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.) @@ -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 @@ -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) @@ -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 @@ -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 "//&