From d7b8e175067dbbd09bdff6e799868a0f056e3d7e Mon Sep 17 00:00:00 2001 From: Marshall Ward Date: Mon, 5 Oct 2020 14:50:38 -0400 Subject: [PATCH 1/3] Dimension: Diagnostic fixes This patch fixes the dimensional scaling of three diagnostics: - grid_Re_Ah - grid_Re_Kh - boundary_forcing_h_tendency The fix to `grid_Re_Ah` and `grid_Re_Kh` removes the `Kh_min` and `Ah_min` parameters, which could not be dimensionally scaled and replaces them with variables in the control structure. The minimum Laplacian and biharmonic viscosites are chosed to estimate the smallest viscosity which is capable to altering the velocity. We use `spacing(1.) * L^2N / T` where `spacing()` is the smallest difference in floating point. We use the minimum harmonic grid spacing, precomputed during initialization, rather than the local value, in order to avoid additional computation and memory access during timestepping. This will changes the value of grid_Re_[AK]h in certain isolated examples of negligible viscosity, but will otherwise avoid the possibility of divide-by-zero while also producing a realizable large value at near-inviscid upper limit. The second fix to `boundary_focing_h_tendency` had its T units set but not its H units. This patch fixes the dimensionality. These are required to fix the dimensional verification errors in PR 1210. --- .../lateral/MOM_hor_visc.F90 | 41 ++++++++++++++----- .../vertical/MOM_diabatic_driver.F90 | 4 +- 2 files changed, 33 insertions(+), 12 deletions(-) diff --git a/src/parameterizations/lateral/MOM_hor_visc.F90 b/src/parameterizations/lateral/MOM_hor_visc.F90 index 76fa656942..a4e2b32b80 100644 --- a/src/parameterizations/lateral/MOM_hor_visc.F90 +++ b/src/parameterizations/lateral/MOM_hor_visc.F90 @@ -93,6 +93,10 @@ module MOM_hor_visc !! depth is shallower than GME_H0 [Z ~> m] real :: GME_efficiency !< The nondimensional prefactor multiplying the GME coefficient [nondim] real :: GME_limiter !< The absolute maximum value the GME coefficient is allowed to take [L2 T-1 ~> m2 s-1]. + real :: min_grid_Kh !< Minimum horizontal Laplacian viscosity used to + !! limit the grid Reynolds number [L2 T-1 ~> m2 s-1] + real :: min_grid_Ah !< Minimun horizontal biharmonic viscosity used to + !! limit grid Reynolds number [L4 T-1 ~> m4 s-1] real ALLOCABLE_, dimension(NIMEM_,NJMEM_) :: Kh_bg_xx !< The background Laplacian viscosity at h points [L2 T-1 ~> m2 s-1]. @@ -358,10 +362,6 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, ! calculation gives the same value as if f were 0 [nondim]. real :: H0_GME ! Depth used to scale down GME coefficient in shallow areas [Z ~> m] real :: KE ! Local kinetic energy [L2 T-2 ~> m2 s-2] - real, parameter :: KH_min = 1.E-30 ! This is the minimun horizontal Laplacian viscosity used to estimate the - ! grid Raynolds number [L2 T-1 ~> m2 s-1] - real, parameter :: AH_min = 1.E-30 ! This is the minimun horizontal Biharmonic viscosity used to estimate the - ! grid Raynolds number [L4 T-1 ~> m4 s-1] logical :: rescale_Kh, legacy_bound logical :: find_FrictWork @@ -877,7 +877,8 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, if (CS%id_grid_Re_Kh>0) then KE = 0.125*((u(I,j,k)+u(I-1,j,k))**2 + (v(i,J,k)+v(i,J-1,k))**2) - grid_Re_Kh(i,j,k) = (sqrt(KE) * sqrt(CS%grid_sp_h2(i,j)))/MAX(Kh,KH_min) + grid_Re_Kh(i,j,k) = (sqrt(KE) * sqrt(CS%grid_sp_h2(i,j))) & + / max(Kh, CS%min_grid_Kh) endif if (CS%id_div_xx_h>0) div_xx_h(i,j,k) = div_xx(i,j) @@ -930,8 +931,9 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, if ((CS%id_Ah_h>0) .or. find_FrictWork .or. CS%debug) Ah_h(i,j,k) = Ah if (CS%id_grid_Re_Ah>0) then - KE = 0.125*((u(I,j,k)+u(I-1,j,k))**2 + (v(i,J,k)+v(i,J-1,k))**2) - grid_Re_Ah(i,j,k) = (sqrt(KE) * CS%grid_sp_h3(i,j))/MAX(Ah,AH_min) + KE = 0.125*((u(I,j,k)+u(I-1,j,k))**2 + (v(i,J,k)+v(i,J-1,k))**2) + grid_Re_Ah(i,j,k) = (sqrt(KE) * CS%grid_sp_h3(i,j)) & + / max(Ah, CS%min_grid_Ah) endif str_xx(i,j) = str_xx(i,j) + Ah * & @@ -1430,6 +1432,8 @@ subroutine hor_visc_init(Time, G, US, param_file, diag, CS, MEKE, ADp) real :: grid_sp_h3 ! Harmonic mean of the squares of the grid^(3/2) [L3 ~> m3] real :: grid_sp_q2 ! spacings at h and q points [L2 ~> m2] real :: grid_sp_q3 ! spacings at h and q points^(3/2) [L3 ~> m3] + real :: min_grid_sp_h2 ! Minimum value of grid_sp_h2 [L2 ~> m2] + real :: min_grid_sp_h4 ! Minimum value of grid_sp_h2**2 [L4 ~> m4] real :: Kh_Limit ! A coefficient [T-1 ~> s-1] used, along with the ! grid spacing, to limit Laplacian viscosity. real :: fmax ! maximum absolute value of f at the four @@ -1718,10 +1722,11 @@ subroutine hor_visc_init(Time, G, US, param_file, diag, CS, MEKE, ADp) "The absolute maximum value the GME coefficient is allowed to take.", & units="m2 s-1", scale=US%m_to_L**2*US%T_to_s, default=1.0e7) endif - if (CS%bound_Kh .or. CS%bound_Ah .or. CS%better_bound_Kh .or. CS%better_bound_Ah) & + if (CS%Laplacian .or. CS%biharmonic) & call get_param(param_file, mdl, "DT", dt, & "The (baroclinic) dynamics time step.", units="s", scale=US%s_to_T, & fail_if_missing=.true.) + Idt = 1.0 / dt if (CS%no_slip .and. CS%biharmonic) & call MOM_error(FATAL,"ERROR: NOSLIP and BIHARMONIC cannot be defined "// & "at the same time in MOM.") @@ -1860,6 +1865,8 @@ subroutine hor_visc_init(Time, G, US, param_file, diag, CS, MEKE, ADp) ! this to be less than 1/3, rather than 1/2 as before. if (CS%bound_Kh .or. CS%bound_Ah) Kh_Limit = 0.3 / (dt*4.0) ! Calculate and store the background viscosity at h-points + + min_grid_sp_h2 = huge(1.) do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 ! Static factors in the Smagorinsky and Leith schemes grid_sp_h2 = (2.0*CS%dx2h(i,j)*CS%dy2h(i,j)) / (CS%dx2h(i,j) + CS%dy2h(i,j)) @@ -1881,7 +1888,9 @@ subroutine hor_visc_init(Time, G, US, param_file, diag, CS, MEKE, ADp) CS%Kh_Max_xx(i,j) = Kh_Limit * grid_sp_h2 CS%Kh_bg_xx(i,j) = MIN(CS%Kh_bg_xx(i,j), CS%Kh_Max_xx(i,j)) endif + min_grid_sp_h2 = min(grid_sp_h2, min_grid_sp_h2) enddo ; enddo + ! Calculate and store the background viscosity at q-points do J=js-1,Jeq ; do I=is-1,Ieq ! Static factors in the Smagorinsky and Leith schemes @@ -1925,6 +1934,8 @@ subroutine hor_visc_init(Time, G, US, param_file, diag, CS, MEKE, ADp) if (CS%better_bound_Ah .or. CS%bound_Ah) Ah_Limit = 0.3 / (dt*64.0) if (CS%Smagorinsky_Ah .and. CS%bound_Coriolis) & BoundCorConst = 1.0 / (5.0*(bound_Cor_vel*bound_Cor_vel)) + + min_grid_sp_h4 = huge(1.) do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 grid_sp_h2 = (2.0*CS%dx2h(i,j)*CS%dy2h(i,j)) / (CS%dx2h(i,j)+CS%dy2h(i,j)) grid_sp_h3 = grid_sp_h2*sqrt(grid_sp_h2) @@ -1949,7 +1960,9 @@ subroutine hor_visc_init(Time, G, US, param_file, diag, CS, MEKE, ADp) CS%Ah_Max_xx(i,j) = Ah_Limit * (grid_sp_h2 * grid_sp_h2) CS%Ah_bg_xx(i,j) = MIN(CS%Ah_bg_xx(i,j), CS%Ah_Max_xx(i,j)) endif + min_grid_sp_h4 = min(grid_sp_h2**2, min_grid_sp_h4) enddo ; enddo + do J=js-1,Jeq ; do I=is-1,Ieq grid_sp_q2 = (2.0*CS%dx2q(I,J)*CS%dy2q(I,J)) / (CS%dx2q(I,J)+CS%dy2q(I,J)) grid_sp_q3 = grid_sp_q2*sqrt(grid_sp_q2) @@ -1975,7 +1988,6 @@ subroutine hor_visc_init(Time, G, US, param_file, diag, CS, MEKE, ADp) endif ! The Laplacian bounds should avoid overshoots when CS%bound_coef < 1. if (CS%Laplacian .and. CS%better_bound_Kh) then - Idt = 1.0 / dt do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 denom = max( & (CS%dy2h(i,j) * CS%DY_dxT(i,j) * (G%IdyCu(I,j) + G%IdyCu(I-1,j)) * & @@ -2004,7 +2016,6 @@ subroutine hor_visc_init(Time, G, US, param_file, diag, CS, MEKE, ADp) ! The biharmonic bounds should avoid overshoots when CS%bound_coef < 0.5, but ! empirically work for CS%bound_coef <~ 1.0 if (CS%biharmonic .and. CS%better_bound_Ah) then - Idt = 1.0 / dt do j=js-1,Jeq+1 ; do I=Isq-1,Ieq+1 u0u(I,j) = (CS%Idxdy2u(I,j)*(CS%dy2h(i+1,j)*CS%DY_dxT(i+1,j)*(G%IdyCu(I+1,j) + G%IdyCu(I,j)) + & CS%dy2h(i,j) * CS%DY_dxT(i,j) * (G%IdyCu(I,j) + G%IdyCu(I-1,j)) ) + & @@ -2104,6 +2115,11 @@ subroutine hor_visc_init(Time, G, US, param_file, diag, CS, MEKE, ADp) 'Biharmonic Horizontal Viscosity at q Points', 'm4 s-1', conversion=US%L_to_m**4*US%s_to_T) CS%id_grid_Re_Ah = register_diag_field('ocean_model', 'grid_Re_Ah', diag%axesTL, Time, & 'Grid Reynolds number for the Biharmonic horizontal viscosity at h points', 'nondim') + + if (CS%id_grid_Re_Ah > 0) & + ! Compute the smallest biharmonic viscosity capable of modifying the + ! velocity at floating point precision. + CS%min_grid_Ah = spacing(1.) * min_grid_sp_h4 * Idt endif if (CS%Laplacian) then CS%id_Kh_h = register_diag_field('ocean_model', 'Khh', diag%axesTL, Time, & @@ -2123,6 +2139,11 @@ subroutine hor_visc_init(Time, G, US, param_file, diag, CS, MEKE, ADp) 'Shearing strain at q Points', 's-1', conversion=US%s_to_T) CS%id_sh_xx_h = register_diag_field('ocean_model', 'sh_xx_h', diag%axesTL, Time, & 'Horizontal tension at h Points', 's-1', conversion=US%s_to_T) + + if (CS%id_grid_Re_Kh > 0) & + ! Compute a smallest Laplacian viscosity capable of modifying the + ! velocity at floating point precision. + CS%min_grid_Kh = spacing(1.) * min_grid_sp_h2 * Idt endif if (CS%use_GME) then CS%id_GME_coeff_h = register_diag_field('ocean_model', 'GME_coeff_h', diag%axesTL, Time, & diff --git a/src/parameterizations/vertical/MOM_diabatic_driver.F90 b/src/parameterizations/vertical/MOM_diabatic_driver.F90 index 62bfeb726d..3d3970f5b0 100644 --- a/src/parameterizations/vertical/MOM_diabatic_driver.F90 +++ b/src/parameterizations/vertical/MOM_diabatic_driver.F90 @@ -3575,8 +3575,8 @@ subroutine diabatic_driver_init(Time, G, GV, US, param_file, useALEalgorithm, di units='m', conversion=GV%H_to_m, v_extensive=.true.) CS%id_boundary_forcing_h_tendency = register_diag_field('ocean_model', & 'boundary_forcing_h_tendency', diag%axesTL, Time, & - 'Cell thickness tendency due to boundary forcing', 'm s-1', conversion=US%s_to_T, & - v_extensive = .true.) + 'Cell thickness tendency due to boundary forcing', 'm s-1', & + conversion=GV%H_to_m*US%s_to_T, v_extensive=.true.) if (CS%id_boundary_forcing_h_tendency > 0) then CS%boundary_forcing_tendency_diag = .true. endif From cb82761df2c54fa1d344e1f66a39fd6dc793b45d Mon Sep 17 00:00:00 2001 From: Marshall Ward Date: Mon, 5 Oct 2020 15:04:17 -0400 Subject: [PATCH 2/3] Bugfix: DT parsing in MOM_hor_visc The single-line if-block for reading the DT parameter did not contain the Idt = 1. / dt update. This patch amends this bug. --- src/parameterizations/lateral/MOM_hor_visc.F90 | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/parameterizations/lateral/MOM_hor_visc.F90 b/src/parameterizations/lateral/MOM_hor_visc.F90 index a4e2b32b80..08b1fec20c 100644 --- a/src/parameterizations/lateral/MOM_hor_visc.F90 +++ b/src/parameterizations/lateral/MOM_hor_visc.F90 @@ -1722,11 +1722,12 @@ subroutine hor_visc_init(Time, G, US, param_file, diag, CS, MEKE, ADp) "The absolute maximum value the GME coefficient is allowed to take.", & units="m2 s-1", scale=US%m_to_L**2*US%T_to_s, default=1.0e7) endif - if (CS%Laplacian .or. CS%biharmonic) & + if (CS%Laplacian .or. CS%biharmonic) then call get_param(param_file, mdl, "DT", dt, & "The (baroclinic) dynamics time step.", units="s", scale=US%s_to_T, & fail_if_missing=.true.) Idt = 1.0 / dt + endif if (CS%no_slip .and. CS%biharmonic) & call MOM_error(FATAL,"ERROR: NOSLIP and BIHARMONIC cannot be defined "// & "at the same time in MOM.") From cbe88ae95796fed9bdf80eeb2f5c903130a584c5 Mon Sep 17 00:00:00 2001 From: Marshall Ward Date: Tue, 6 Oct 2020 14:49:31 -0400 Subject: [PATCH 3/3] min_grid_sp_h[24] across PEs An additional `min_across_PEs` call is added to the min_grid_sp_h[24] calculations, to accommodate for parallel builds and ensure consistency across layouts. --- src/parameterizations/lateral/MOM_hor_visc.F90 | 3 +++ 1 file changed, 3 insertions(+) diff --git a/src/parameterizations/lateral/MOM_hor_visc.F90 b/src/parameterizations/lateral/MOM_hor_visc.F90 index 08b1fec20c..3c63564a30 100644 --- a/src/parameterizations/lateral/MOM_hor_visc.F90 +++ b/src/parameterizations/lateral/MOM_hor_visc.F90 @@ -4,6 +4,7 @@ module MOM_hor_visc ! This file is part of MOM6. See LICENSE.md for the license. use MOM_checksums, only : hchksum, Bchksum +use MOM_coms, only : min_across_PEs use MOM_diag_mediator, only : post_data, register_diag_field, safe_alloc_ptr use MOM_diag_mediator, only : diag_ctrl, time_type use MOM_domains, only : pass_var, CORNER, pass_vector, AGRID, BGRID_NE @@ -1891,6 +1892,7 @@ subroutine hor_visc_init(Time, G, US, param_file, diag, CS, MEKE, ADp) endif min_grid_sp_h2 = min(grid_sp_h2, min_grid_sp_h2) enddo ; enddo + call min_across_PEs(min_grid_sp_h2) ! Calculate and store the background viscosity at q-points do J=js-1,Jeq ; do I=is-1,Ieq @@ -1963,6 +1965,7 @@ subroutine hor_visc_init(Time, G, US, param_file, diag, CS, MEKE, ADp) endif min_grid_sp_h4 = min(grid_sp_h2**2, min_grid_sp_h4) enddo ; enddo + call min_across_PEs(min_grid_sp_h4) do J=js-1,Jeq ; do I=is-1,Ieq grid_sp_q2 = (2.0*CS%dx2q(I,J)*CS%dy2q(I,J)) / (CS%dx2q(I,J)+CS%dy2q(I,J))