From 184d24f6055bb4e2c4470b92618295cbd4eed85c Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Tue, 26 Oct 2021 04:26:21 -0400 Subject: [PATCH 1/2] Eliminate GET_ALL_PARAMS in hor_visc_init Added do_not_log arguments to get_param calls in MOM_hor_visc.F90 that are only used conditionally, and eliminated the unlogged GET_ALL_PARAMS runtime parameter and get_all variable in hor_visc_init(). By design, all logging of parameters after this commit is identical to before, even for variables that are inactive and therefore should not be logged. In several places, there were some problems, mostly with the GME code, that have been noted in comments marked with '###'. Also cleaned up the code alignment and eliminated unneeded temporary variables in a few places in hor_visc(). All solutions are bitwise identical, and no output is changed. --- .../lateral/MOM_hor_visc.F90 | 339 +++++++++--------- 1 file changed, 170 insertions(+), 169 deletions(-) diff --git a/src/parameterizations/lateral/MOM_hor_visc.F90 b/src/parameterizations/lateral/MOM_hor_visc.F90 index 15e8415474..35a1f511b3 100644 --- a/src/parameterizations/lateral/MOM_hor_visc.F90 +++ b/src/parameterizations/lateral/MOM_hor_visc.F90 @@ -338,8 +338,6 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, ! viscosity. Here set equal to nondimensional Laplacian Leith constant. ! This is set equal to zero if modified Leith is not used. real :: Shear_mag_bc ! Shear_mag value in backscatter [T-1 ~> s-1] - real :: sh_xx_sq ! Square of tension (sh_xx) [T-2 ~> s-2] - real :: sh_xy_sq ! Square of shearing strain (sh_xy) [T-2 ~> s-2] real :: h2uq, h2vq ! temporary variables [H2 ~> m2 or kg2 m-4]. real :: hu, hv ! Thicknesses interpolated by arithmetic means to corner ! points; these are first interpolated to u or v velocity @@ -573,7 +571,7 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, !$OMP grad_vel_mag_h, grad_vel_mag_q, & !$OMP grad_vel_mag_bt_h, grad_vel_mag_bt_q, grad_d2vel_mag_h, & !$OMP meke_res_fn, Shear_mag, Shear_mag_bc, vert_vort_mag, h_min, hrat_min, visc_bound_rem, & - !$OMP sh_xx_sq, sh_xy_sq, grid_Ah, grid_Kh, d_Del2u, d_Del2v, d_str, & + !$OMP grid_Ah, grid_Kh, d_Del2u, d_Del2v, d_str, & !$OMP Kh, Ah, AhSm, AhLth, local_strain, Sh_F_pow, & !$OMP dDel2vdx, dDel2udy, DY_dxCv, DX_dyCu, Del2vort_q, Del2vort_h, KE, & !$OMP h2uq, h2vq, hu, hv, hq, FatH, RoScl, GME_coeff & @@ -819,7 +817,7 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, DX_dyBu = G%dxBu(I,J) * G%IdyBu(I,J) Del2vort_q(I,J) = DY_dxBu * (vort_xy_dx(i+1,J) * G%IdyCv(i+1,J) - vort_xy_dx(i,J) * G%IdyCv(i,J)) + & - DX_dyBu * (vort_xy_dy(I,j+1) * G%IdyCu(I,j+1) - vort_xy_dy(I,j) * G%IdyCu(I,j)) + DX_dyBu * (vort_xy_dy(I,j+1) * G%IdyCu(I,j+1) - vort_xy_dy(I,j) * G%IdyCu(I,j)) enddo ; enddo do J=Jsq,Jeq+1 ; do I=Isq,Ieq+1 Del2vort_h(i,j) = 0.25*(Del2vort_q(I,J) + Del2vort_q(I-1,J) + Del2vort_q(I,J-1) + Del2vort_q(I-1,J-1)) @@ -837,12 +835,12 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, ! Magnitude of divergence gradient do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 - grad_div_mag_h(i,j) =sqrt((0.5*(div_xx_dx(I,j) + div_xx_dx(I-1,j)))**2 + & - (0.5 * (div_xx_dy(i,J) + div_xx_dy(i,J-1)))**2) + grad_div_mag_h(i,j) = sqrt((0.5*(div_xx_dx(I,j) + div_xx_dx(I-1,j)))**2 + & + (0.5*(div_xx_dy(i,J) + div_xx_dy(i,J-1)))**2) enddo ; enddo do j=Jsq-1,Jeq+1 ; do i=Isq-1,Ieq+1 - grad_div_mag_q(I,J) =sqrt((0.5*(div_xx_dx(I,j) + div_xx_dx(I,j+1)))**2 + & - (0.5 * (div_xx_dy(i,J) + div_xx_dy(i+1,J)))**2) + grad_div_mag_q(I,J) = sqrt((0.5*(div_xx_dx(I,j) + div_xx_dx(I,j+1)))**2 + & + (0.5*(div_xx_dy(i,J) + div_xx_dy(i+1,J)))**2) enddo ; enddo else @@ -902,12 +900,9 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, if ((CS%Smagorinsky_Kh) .or. (CS%Smagorinsky_Ah)) then do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 - sh_xx_sq = sh_xx(i,j) * sh_xx(i,j) - sh_xy_sq = 0.25 * ( & - (sh_xy(I-1,J-1) * sh_xy(I-1,J-1) + sh_xy(I,J) * sh_xy(I,J)) & - + (sh_xy(I-1,J) * sh_xy(I-1,J) + sh_xy(I,J-1) * sh_xy(I,J-1)) & - ) - Shear_mag(i,j) = sqrt(sh_xx_sq + sh_xy_sq) + Shear_mag(i,j) = sqrt( sh_xx(i,j)**2 + & + 0.25*((sh_xy(I-1,J-1)**2 + sh_xy(I,J)**2) + & + (sh_xy(I-1,J)**2 + sh_xy(I,J-1)**2)) ) enddo ; enddo endif @@ -1184,12 +1179,9 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, if ((CS%Smagorinsky_Kh) .or. (CS%Smagorinsky_Ah)) then do J=js-1,Jeq ; do I=is-1,Ieq - sh_xy_sq = sh_xy(I,J) * sh_xy(I,J) - sh_xx_sq = 0.25 * ( & - (sh_xx(i,j) * sh_xx(i,j) + sh_xx(i+1,j+1) * sh_xx(i+1,j+1)) & - + (sh_xx(i,j+1) * sh_xx(i,j+1) + sh_xx(i+1,j) * sh_xx(i+1,j)) & - ) - Shear_mag(i,j) = sqrt(sh_xy_sq + sh_xx_sq) + Shear_mag(I,J) = sqrt( sh_xy(I,J)**2 + & + 0.25*((sh_xx(i,j)**2 + sh_xx(i+1,j+1)**2) + & + (sh_xx(i,j+1)**2 + sh_xx(i+1,j)**2)) ) enddo ; enddo endif @@ -1434,6 +1426,9 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, endif if (CS%use_GME) then + !### This call to get the 3-d GME diffusivity arrays and the subsequent blocking halo update + ! should occur outside of the k-loop, and perhaps the halo update should occur outside of + ! this routine altogether! call thickness_diffuse_get_KH(TD, KH_u_GME, KH_v_GME, G, GV) call pass_vector(KH_u_GME, KH_v_GME, G%Domain) @@ -1458,6 +1453,9 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, enddo ; enddo ! Applying GME diagonal term. This is linear and the arguments can be rescaled. + !### This smoothing is only applied at computational grid points, but is used in extra halo points! + !### There are blocking halo updates in the smooth_GME routines, which could be avoided by expanding + ! the loop ranges by a point in the code setting str_xx_GME and str_xy_GME a few lines above. call smooth_GME(CS, G, GME_flux_h=str_xx_GME) call smooth_GME(CS, G, GME_flux_q=str_xy_GME) @@ -1809,9 +1807,6 @@ subroutine hor_visc_init(Time, G, GV, US, param_file, diag, CS, MEKE, ADp) real :: Kh_sin_lat ! Amplitude of latitudinally dependent viscosity [L2 T-1 ~> m2 s-1] real :: Kh_pwr_of_sine ! Power used to raise sin(lat) when using Kh_sin_lat logical :: bound_Cor_def ! parameter setting of BOUND_CORIOLIS - logical :: get_all ! If true, read and log all parameters, regardless of - ! whether they are used, to enable spell-checking of - ! valid parameters. logical :: split ! If true, use the split time stepping scheme. ! If false and USE_GME = True, issue a FATAL error. logical :: use_MEKE ! If true, the MEKE parameterization is in use. @@ -1824,8 +1819,8 @@ subroutine hor_visc_init(Time, G, GV, US, param_file, diag, CS, MEKE, ADp) integer :: is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz integer :: isd, ied, jsd, jed, IsdB, IedB, JsdB, JedB integer :: i, j -! This include declares and sets the variable "version". -#include "version_variable.h" + ! This include declares and sets the variable "version". +# include "version_variable.h" character(len=40) :: mdl = "MOM_hor_visc" ! module name is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke Isq = G%IscB ; Ieq = G%IecB ; Jsq = G%JscB ; Jeq = G%JecB @@ -1840,22 +1835,8 @@ subroutine hor_visc_init(Time, G, GV, US, param_file, diag, CS, MEKE, ADp) CS%diag => diag ! Read parameters and write them to the model log. call log_version(param_file, mdl, version, "") - ! It is not clear whether all of these initialization lines are needed for the - ! cases where the corresponding parameters are not read. - CS%bound_Kh = .false. ; CS%better_bound_Kh = .false. ; CS%Smagorinsky_Kh = .false. ; CS%Leith_Kh = .false. - CS%bound_Ah = .false. ; CS%better_bound_Ah = .false. ; CS%Smagorinsky_Ah = .false. ; CS%Leith_Ah = .false. - CS%use_QG_Leith_visc = .false. - CS%bound_Coriolis = .false. - CS%Modified_Leith = .false. - CS%dynamic_aniso = .false. - Kh = 0.0 ; Ah = 0.0 - ! These initialization lines are needed because they are used even in cases where they are not read. - CS%anisotropic = .false. - CS%res_scale_MEKE = .false. - - ! If GET_ALL_PARAMS is true, all parameters are read in all cases to enable - ! parameter spelling checks. - call get_param(param_file, mdl, "GET_ALL_PARAMS", get_all, default=.false.) + + ! All parameters are read in all cases to enable parameter spelling checks. call get_param(param_file, mdl, "DEFAULT_2018_ANSWERS", default_2018_answers, & "This sets the default value for the various _2018_ANSWERS parameters.", & default=.false.) @@ -1867,185 +1848,202 @@ subroutine hor_visc_init(Time, G, GV, US, param_file, diag, CS, MEKE, ADp) call get_param(param_file, mdl, "LAPLACIAN", CS%Laplacian, & "If true, use a Laplacian horizontal viscosity.", & default=.false.) - if (CS%Laplacian .or. get_all) then - call get_param(param_file, mdl, "KH", Kh, & + + call get_param(param_file, mdl, "KH", Kh, & "The background Laplacian horizontal viscosity.", & - units = "m2 s-1", default=0.0, scale=US%m_to_L**2*US%T_to_s) - call get_param(param_file, mdl, "KH_BG_MIN", CS%Kh_bg_min, & + units="m2 s-1", default=0.0, scale=US%m_to_L**2*US%T_to_s, & + do_not_log=.not.CS%Laplacian) + call get_param(param_file, mdl, "KH_BG_MIN", CS%Kh_bg_min, & "The minimum value allowed for Laplacian horizontal viscosity, KH.", & - units = "m2 s-1", default=0.0, scale=US%m_to_L**2*US%T_to_s) - call get_param(param_file, mdl, "KH_VEL_SCALE", Kh_vel_scale, & + units="m2 s-1", default=0.0, scale=US%m_to_L**2*US%T_to_s, & + do_not_log=.not.CS%Laplacian) + call get_param(param_file, mdl, "KH_VEL_SCALE", Kh_vel_scale, & "The velocity scale which is multiplied by the grid "//& "spacing to calculate the Laplacian viscosity. "//& "The final viscosity is the largest of this scaled "//& "viscosity, the Smagorinsky and Leith viscosities, and KH.", & - units="m s-1", default=0.0, scale=US%m_s_to_L_T) - call get_param(param_file, mdl, "KH_SIN_LAT", Kh_sin_lat, & + units="m s-1", default=0.0, scale=US%m_s_to_L_T, & + do_not_log=.not.CS%Laplacian) + call get_param(param_file, mdl, "KH_SIN_LAT", Kh_sin_lat, & "The amplitude of a latitudinally-dependent background "//& "viscosity of the form KH_SIN_LAT*(SIN(LAT)**KH_PWR_OF_SINE).", & - units = "m2 s-1", default=0.0, scale=US%m_to_L**2*US%T_to_s) - if (Kh_sin_lat>0. .or. get_all) & - call get_param(param_file, mdl, "KH_PWR_OF_SINE", Kh_pwr_of_sine, & + units="m2 s-1", default=0.0, scale=US%m_to_L**2*US%T_to_s, & + do_not_log=.not.CS%Laplacian) + call get_param(param_file, mdl, "KH_PWR_OF_SINE", Kh_pwr_of_sine, & "The power used to raise SIN(LAT) when using a latitudinally "//& "dependent background viscosity.", & - units = "nondim", default=4.0) - call get_param(param_file, mdl, "SMAGORINSKY_KH", CS%Smagorinsky_Kh, & + units="nondim", default=4.0, & + do_not_log=.not.(CS%Laplacian .and. (Kh_sin_lat>0.)) ) + call get_param(param_file, mdl, "SMAGORINSKY_KH", CS%Smagorinsky_Kh, & "If true, use a Smagorinsky nonlinear eddy viscosity.", & - default=.false.) - if (CS%Smagorinsky_Kh .or. get_all) & - call get_param(param_file, mdl, "SMAG_LAP_CONST", Smag_Lap_const, & + default=.false., do_not_log=.not.CS%Laplacian) + if (.not.CS%Laplacian) CS%Smagorinsky_Kh = .false. + call get_param(param_file, mdl, "SMAG_LAP_CONST", Smag_Lap_const, & "The nondimensional Laplacian Smagorinsky constant, "//& "often 0.15.", units="nondim", default=0.0, & - fail_if_missing = CS%Smagorinsky_Kh) - call get_param(param_file, mdl, "LEITH_KH", CS%Leith_Kh, & + fail_if_missing=CS%Smagorinsky_Kh, do_not_log=.not.CS%Smagorinsky_Kh) + call get_param(param_file, mdl, "LEITH_KH", CS%Leith_Kh, & "If true, use a Leith nonlinear eddy viscosity.", & - default=.false.) - ! This call duplicates one that occurs 26 lines later, and is probably unneccessary. - call get_param(param_file, mdl, "MODIFIED_LEITH", CS%Modified_Leith, & + default=.false., do_not_log=.not.CS%Laplacian) + if (.not.CS%Laplacian) CS%Leith_Kh = .false. + ! This call duplicates one that occurs 26 lines later, and is probably unneccessary. + call get_param(param_file, mdl, "MODIFIED_LEITH", CS%Modified_Leith, & "If true, add a term to Leith viscosity which is "//& "proportional to the gradient of divergence.", & - default=.false.) - call get_param(param_file, mdl, "USE_MEKE", use_MEKE, & + default=.false., do_not_log=.not.CS%Laplacian) !### (.not.CS%Leith_Kh)? + call get_param(param_file, mdl, "USE_MEKE", use_MEKE, & default=.false., do_not_log=.true.) - call get_param(param_file, mdl, "RES_SCALE_MEKE_VISC", CS%res_scale_MEKE, & + call get_param(param_file, mdl, "RES_SCALE_MEKE_VISC", CS%res_scale_MEKE, & "If true, the viscosity contribution from MEKE is scaled by "//& - "the resolution function.", default=.false., do_not_log=.not.use_MEKE) - if (.not.use_MEKE) CS%res_scale_MEKE = .false. - if (CS%Leith_Kh .or. get_all) then - call get_param(param_file, mdl, "LEITH_LAP_CONST", Leith_Lap_const, & + "the resolution function.", default=.false., & + do_not_log=.not.(CS%Laplacian.and.use_MEKE)) + if (.not.(CS%Laplacian.and.use_MEKE)) CS%res_scale_MEKE = .false. + + call get_param(param_file, mdl, "LEITH_LAP_CONST", Leith_Lap_const, & "The nondimensional Laplacian Leith constant, "//& "often set to 1.0", units="nondim", default=0.0, & - fail_if_missing = CS%Leith_Kh) - call get_param(param_file, mdl, "USE_QG_LEITH_VISC", CS%use_QG_Leith_visc, & + fail_if_missing=CS%Leith_Kh, do_not_log=.not.CS%Leith_Kh) + call get_param(param_file, mdl, "USE_QG_LEITH_VISC", CS%use_QG_Leith_visc, & "If true, use QG Leith nonlinear eddy viscosity.", & - default=.false.) - if (CS%use_QG_Leith_visc .and. .not. CS%Leith_Kh) call MOM_error(FATAL, & + default=.false., do_not_log=.not.CS%Leith_Kh) + if (CS%use_QG_Leith_visc .and. .not. CS%Leith_Kh) call MOM_error(FATAL, & "MOM_hor_visc.F90, hor_visc_init:"//& "LEITH_KH must be True when USE_QG_LEITH_VISC=True.") - endif - if (CS%Leith_Kh .or. CS%Leith_Ah .or. get_all) then - call get_param(param_file, mdl, "USE_BETA_IN_LEITH", CS%use_beta_in_Leith, & + + !### The following two get_param_calls need to occur after Leith_Ah is read, but for now it replciates prior code. + CS%Leith_Ah = .false. + call get_param(param_file, mdl, "USE_BETA_IN_LEITH", CS%use_beta_in_Leith, & "If true, include the beta term in the Leith nonlinear eddy viscosity.", & - default=CS%Leith_Kh) - call get_param(param_file, mdl, "MODIFIED_LEITH", CS%modified_Leith, & + default=CS%Leith_Kh, do_not_log=.not.(CS%Leith_Kh .or. CS%Leith_Ah) ) + call get_param(param_file, mdl, "MODIFIED_LEITH", CS%modified_Leith, & "If true, add a term to Leith viscosity which is "//& "proportional to the gradient of divergence.", & - default=.false.) - endif - call get_param(param_file, mdl, "BOUND_KH", CS%bound_Kh, & + default=.false., do_not_log=.not.(CS%Leith_Kh .or. CS%Leith_Ah) ) + + call get_param(param_file, mdl, "BOUND_KH", CS%bound_Kh, & "If true, the Laplacian coefficient is locally limited "//& - "to be stable.", default=.true.) - call get_param(param_file, mdl, "BETTER_BOUND_KH", CS%better_bound_Kh, & + "to be stable.", default=.true., do_not_log=.not.CS%Laplacian) + call get_param(param_file, mdl, "BETTER_BOUND_KH", CS%better_bound_Kh, & "If true, the Laplacian coefficient is locally limited "//& "to be stable with a better bounding than just BOUND_KH.", & - default=CS%bound_Kh) - call get_param(param_file, mdl, "ANISOTROPIC_VISCOSITY", CS%anisotropic, & + default=CS%bound_Kh, do_not_log=.not.CS%Laplacian) + if (.not.CS%Laplacian) CS%bound_Kh = .false. + if (.not.CS%Laplacian) CS%better_bound_Kh = .false. + call get_param(param_file, mdl, "ANISOTROPIC_VISCOSITY", CS%anisotropic, & "If true, allow anistropic viscosity in the Laplacian "//& - "horizontal viscosity.", default=.false.) - call get_param(param_file, mdl, "ADD_LES_VISCOSITY", CS%add_LES_viscosity, & + "horizontal viscosity.", default=.false., & + do_not_log=.not.CS%Laplacian) + if (.not.CS%Laplacian) CS%anisotropic = .false. ! This replicates the prior code, but is it intended? + call get_param(param_file, mdl, "ADD_LES_VISCOSITY", CS%add_LES_viscosity, & "If true, adds the viscosity from Smagorinsky and Leith to the "//& - "background viscosity instead of taking the maximum.", default=.false.) - endif - if (CS%anisotropic .or. get_all) then - call get_param(param_file, mdl, "KH_ANISO", CS%Kh_aniso, & + "background viscosity instead of taking the maximum.", default=.false., & + do_not_log=.not.CS%Laplacian) + + call get_param(param_file, mdl, "KH_ANISO", CS%Kh_aniso, & "The background Laplacian anisotropic horizontal viscosity.", & - units = "m2 s-1", default=0.0, scale=US%m_to_L**2*US%T_to_s) - call get_param(param_file, mdl, "ANISOTROPIC_MODE", aniso_mode, & + units="m2 s-1", default=0.0, scale=US%m_to_L**2*US%T_to_s, & + do_not_log=.not.CS%anisotropic) + call get_param(param_file, mdl, "ANISOTROPIC_MODE", aniso_mode, & "Selects the mode for setting the direction of anistropy.\n"//& "\t 0 - Points along the grid i-direction.\n"//& "\t 1 - Points towards East.\n"//& "\t 2 - Points along the flow direction, U/|U|.", & - default=0) - select case (aniso_mode) - case (0) - call get_param(param_file, mdl, "ANISO_GRID_DIR", aniso_grid_dir, & - "The vector pointing in the direction of anistropy for "//& - "horizont viscosity. n1,n2 are the i,j components relative "//& - "to the grid.", units = "nondim", fail_if_missing=.true.) - case (1) - call get_param(param_file, mdl, "ANISO_GRID_DIR", aniso_grid_dir, & - "The vector pointing in the direction of anistropy for "//& - "horizont viscosity. n1,n2 are the i,j components relative "//& - "to the spherical coordinates.", units = "nondim", fail_if_missing=.true.) - end select + default=0, do_not_log=.not.CS%anisotropic) + if (aniso_mode == 0) then + call get_param(param_file, mdl, "ANISO_GRID_DIR", aniso_grid_dir, & + "The vector pointing in the direction of anistropy for horizontal viscosity. "//& + "n1,n2 are the i,j components relative to the grid.", & + units="nondim", fail_if_missing=CS%anisotropic, do_not_log=.not.CS%anisotropic) + elseif (aniso_mode == 1) then + call get_param(param_file, mdl, "ANISO_GRID_DIR", aniso_grid_dir, & + "The vector pointing in the direction of anistropy for horizontal viscosity. "//& + "n1,n2 are the i,j components relative to the spherical coordinates.", & + units="nondim", fail_if_missing=CS%anisotropic, do_not_log=.not.CS%anisotropic) + else + call get_param(param_file, mdl, "ANISO_GRID_DIR", aniso_grid_dir, & + "The vector pointing in the direction of anistropy for horizontal viscosity.", & + units="nondim", fail_if_missing=.false., do_not_log=.true.) endif + call get_param(param_file, mdl, "BIHARMONIC", CS%biharmonic, & "If true, use a biharmonic horizontal viscosity. "//& "BIHARMONIC may be used with LAPLACIAN.", & default=.true.) - if (CS%biharmonic .or. get_all) then - call get_param(param_file, mdl, "AH", Ah, & + call get_param(param_file, mdl, "AH", Ah, & "The background biharmonic horizontal viscosity.", & - units = "m4 s-1", default=0.0, scale=US%m_to_L**4*US%T_to_s) - call get_param(param_file, mdl, "AH_VEL_SCALE", Ah_vel_scale, & + units="m4 s-1", default=0.0, scale=US%m_to_L**4*US%T_to_s, & + do_not_log=.not.CS%biharmonic) + call get_param(param_file, mdl, "AH_VEL_SCALE", Ah_vel_scale, & "The velocity scale which is multiplied by the cube of "//& "the grid spacing to calculate the biharmonic viscosity. "//& "The final viscosity is the largest of this scaled "//& "viscosity, the Smagorinsky and Leith viscosities, and AH.", & - units="m s-1", default=0.0, scale=US%m_s_to_L_T) - call get_param(param_file, mdl, "AH_TIME_SCALE", Ah_time_scale, & + units="m s-1", default=0.0, scale=US%m_s_to_L_T, do_not_log=.not.CS%biharmonic) + call get_param(param_file, mdl, "AH_TIME_SCALE", Ah_time_scale, & "A time scale whose inverse is multiplied by the fourth "//& "power of the grid spacing to calculate biharmonic viscosity. "//& "The final viscosity is the largest of all viscosity "//& "formulations in use. 0.0 means that it's not used.", & - units="s", default=0.0, scale=US%s_to_T) - call get_param(param_file, mdl, "SMAGORINSKY_AH", CS%Smagorinsky_Ah, & + units="s", default=0.0, scale=US%s_to_T, do_not_log=.not.CS%biharmonic) + call get_param(param_file, mdl, "SMAGORINSKY_AH", CS%Smagorinsky_Ah, & "If true, use a biharmonic Smagorinsky nonlinear eddy "//& - "viscosity.", default=.false.) - call get_param(param_file, mdl, "LEITH_AH", CS%Leith_Ah, & + "viscosity.", default=.false., do_not_log=.not.CS%biharmonic) + if (.not.CS%biharmonic) CS%Smagorinsky_Ah = .false. + call get_param(param_file, mdl, "LEITH_AH", CS%Leith_Ah, & "If true, use a biharmonic Leith nonlinear eddy "//& - "viscosity.", default=.false.) - call get_param(param_file, mdl, "BOUND_AH", CS%bound_Ah, & + "viscosity.", default=.false., do_not_log=.not.CS%biharmonic) + if (.not.CS%biharmonic) CS%Leith_Ah = .false. + call get_param(param_file, mdl, "BOUND_AH", CS%bound_Ah, & "If true, the biharmonic coefficient is locally limited "//& - "to be stable.", default=.true.) - call get_param(param_file, mdl, "BETTER_BOUND_AH", CS%better_bound_Ah, & + "to be stable.", default=.true., do_not_log=.not.CS%biharmonic) + call get_param(param_file, mdl, "BETTER_BOUND_AH", CS%better_bound_Ah, & "If true, the biharmonic coefficient is locally limited "//& "to be stable with a better bounding than just BOUND_AH.", & - default=CS%bound_Ah) - call get_param(param_file, mdl, "RE_AH", CS%Re_Ah, & + default=CS%bound_Ah, do_not_log=.not.CS%biharmonic) + if (.not.CS%biharmonic) CS%bound_Ah = .false. + if (.not.CS%biharmonic) CS%better_bound_Ah = .false. + call get_param(param_file, mdl, "RE_AH", CS%Re_Ah, & "If nonzero, the biharmonic coefficient is scaled "//& "so that the biharmonic Reynolds number is equal to this.", & - units="nondim", default=0.0) + units="nondim", default=0.0, do_not_log=.not.CS%biharmonic) - if (CS%Smagorinsky_Ah .or. get_all) then - call get_param(param_file, mdl, "SMAG_BI_CONST",Smag_bi_const, & + call get_param(param_file, mdl, "SMAG_BI_CONST",Smag_bi_const, & "The nondimensional biharmonic Smagorinsky constant, "//& "typically 0.015 - 0.06.", units="nondim", default=0.0, & - fail_if_missing = CS%Smagorinsky_Ah) - call get_param(param_file, mdl, "BOUND_CORIOLIS", bound_Cor_def, default=.false.) - call get_param(param_file, mdl, "BOUND_CORIOLIS_BIHARM", CS%bound_Coriolis, & + fail_if_missing=CS%Smagorinsky_Ah, do_not_log=.not.CS%Smagorinsky_Ah) + + call get_param(param_file, mdl, "BOUND_CORIOLIS", bound_Cor_def, default=.false.) + call get_param(param_file, mdl, "BOUND_CORIOLIS_BIHARM", CS%bound_Coriolis, & "If true use a viscosity that increases with the square "//& "of the velocity shears, so that the resulting viscous "//& "drag is of comparable magnitude to the Coriolis terms "//& "when the velocity differences between adjacent grid "//& "points is 0.5*BOUND_CORIOLIS_VEL. The default is the "//& - "value of BOUND_CORIOLIS (or false).", default=bound_Cor_def) - if (CS%bound_Coriolis .or. get_all) then - call get_param(param_file, mdl, "MAXVEL", maxvel, default=3.0e8) - bound_Cor_vel = maxvel - call get_param(param_file, mdl, "BOUND_CORIOLIS_VEL", bound_Cor_vel, & + "value of BOUND_CORIOLIS (or false).", default=bound_Cor_def, & + do_not_log=.not.CS%Smagorinsky_Ah) + if (.not.CS%Smagorinsky_Ah) CS%bound_Coriolis = .false. + call get_param(param_file, mdl, "MAXVEL", maxvel, default=3.0e8) + call get_param(param_file, mdl, "BOUND_CORIOLIS_VEL", bound_Cor_vel, & "The velocity scale at which BOUND_CORIOLIS_BIHARM causes "//& "the biharmonic drag to have comparable magnitude to the "//& "Coriolis acceleration. The default is set by MAXVEL.", & - units="m s-1", default=maxvel, scale=US%m_s_to_L_T) - endif - endif - if (CS%Leith_Ah .or. get_all) & - call get_param(param_file, mdl, "LEITH_BI_CONST", Leith_bi_const, & + units="m s-1", default=maxvel, scale=US%m_s_to_L_T, & + do_not_log=.not.(CS%Smagorinsky_Ah .and. CS%bound_Coriolis)) + + call get_param(param_file, mdl, "LEITH_BI_CONST", Leith_bi_const, & "The nondimensional biharmonic Leith constant, "//& "typical values are thus far undetermined.", units="nondim", default=0.0, & - fail_if_missing = CS%Leith_Ah) - endif + fail_if_missing=CS%Leith_Ah, do_not_log=.not.CS%Leith_Ah) + call get_param(param_file, mdl, "USE_LAND_MASK_FOR_HVISC", CS%use_land_mask, & "If true, use Use the land mask for the computation of thicknesses "//& "at velocity locations. This eliminates the dependence on arbitrary "//& "values over land or outside of the domain.", default=.true.) - if (CS%better_bound_Ah .or. CS%better_bound_Kh .or. get_all) & - call get_param(param_file, mdl, "HORVISC_BOUND_COEF", CS%bound_coef, & + call get_param(param_file, mdl, "HORVISC_BOUND_COEF", CS%bound_coef, & "The nondimensional coefficient of the ratio of the "//& "viscosity bounds to the theoretical maximum for "//& "stability without considering other terms.", units="nondim", & - default=0.8) + default=0.8, do_not_log=.not.(CS%better_bound_Ah .or. CS%better_bound_Kh)) call get_param(param_file, mdl, "NOSLIP", CS%no_slip, & "If true, no slip boundary conditions are used; otherwise "//& "free slip boundary conditions are assumed. The "//& @@ -2056,34 +2054,31 @@ subroutine hor_visc_init(Time, G, GV, US, param_file, diag, CS, MEKE, ADp) call get_param(param_file, mdl, "USE_KH_BG_2D", CS%use_Kh_bg_2d, & "If true, read a file containing 2-d background harmonic "//& "viscosities. The final viscosity is the maximum of the other "//& - "terms and this background value.", default=.false.) - if (CS%use_Kh_bg_2d) then - call get_param(param_file, mdl, "KH_BG_2D_BUG", CS%Kh_bg_2d_bug, & + "terms and this background value.", default=.false.) ! ###do_not_log=.not.CS%Laplacian? + call get_param(param_file, mdl, "KH_BG_2D_BUG", CS%Kh_bg_2d_bug, & "If true, retain an answer-changing horizontal indexing bug in setting "//& "the corner-point viscosities when USE_KH_BG_2D=True. This is"//& - "not recommended.", default=.false.) - endif + "not recommended.", default=.false., do_not_log=.not.CS%use_Kh_bg_2d) call get_param(param_file, mdl, "USE_GME", CS%use_GME, & "If true, use the GM+E backscatter scheme in association \n"//& "with the Gent and McWilliams parameterization.", default=.false.) - if (CS%use_GME) then - call get_param(param_file, mdl, "SPLIT", split, & - "Use the split time stepping if true.", default=.true., & - do_not_log=.true.) - if (.not. split) call MOM_error(FATAL,"ERROR: Currently, USE_GME = True "// & + call get_param(param_file, mdl, "SPLIT", split, & + "Use the split time stepping if true.", default=.true., do_not_log=.true.) + if (CS%use_GME .and. .not.split) call MOM_error(FATAL,"ERROR: Currently, USE_GME = True "// & "cannot be used with SPLIT=False.") - call get_param(param_file, mdl, "GME_H0", CS%GME_h0, & + call get_param(param_file, mdl, "GME_H0", CS%GME_h0, & "The strength of GME tapers quadratically to zero when the bathymetric "//& - "depth is shallower than GME_H0.", units="m", scale=US%m_to_Z, & - default=1000.0) - call get_param(param_file, mdl, "GME_EFFICIENCY", CS%GME_efficiency, & + "depth is shallower than GME_H0.", & + units="m", scale=US%m_to_Z, default=1000.0, do_not_log=.not.CS%use_GME) + call get_param(param_file, mdl, "GME_EFFICIENCY", CS%GME_efficiency, & "The nondimensional prefactor multiplying the GME coefficient.", & - units="nondim", default=1.0) - call get_param(param_file, mdl, "GME_LIMITER", CS%GME_limiter, & + units="nondim", default=1.0, do_not_log=.not.CS%use_GME) + call get_param(param_file, mdl, "GME_LIMITER", CS%GME_limiter, & "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 + units="m2 s-1", scale=US%m_to_L**2*US%T_to_s, default=1.0e7, & + do_not_log=.not.CS%use_GME) + 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, & @@ -2128,6 +2123,8 @@ subroutine hor_visc_init(Time, G, GV, US, param_file, diag, CS, MEKE, ADp) endif ALLOC_(CS%reduction_xx(isd:ied,jsd:jed)) ; CS%reduction_xx(:,:) = 0.0 ALLOC_(CS%reduction_xy(IsdB:IedB,JsdB:JedB)) ; CS%reduction_xy(:,:) = 0.0 + + CS%dynamic_aniso = .false. if (CS%anisotropic) then ALLOC_(CS%n1n2_h(isd:ied,jsd:jed)) ; CS%n1n2_h(:,:) = 0.0 ALLOC_(CS%n1n1_m_n2n2_h(isd:ied,jsd:jed)) ; CS%n1n1_m_n2n2_h(:,:) = 0.0 @@ -2145,13 +2142,14 @@ subroutine hor_visc_init(Time, G, GV, US, param_file, diag, CS, MEKE, ADp) "Runtime parameter ANISOTROPIC_MODE is out of range.") end select endif - if (CS%use_Kh_bg_2d) then - ALLOC_(CS%Kh_bg_2d(isd:ied,jsd:jed)) ; CS%Kh_bg_2d(:,:) = 0.0 - call get_param(param_file, mdl, "KH_BG_2D_FILENAME", filename, & + + call get_param(param_file, mdl, "KH_BG_2D_FILENAME", filename, & 'The filename containing a 2d map of "Kh".', & - default='KH_background_2d.nc') + default='KH_background_2d.nc', do_not_log=.not.CS%use_Kh_bg_2d) + if (CS%use_Kh_bg_2d) then call get_param(param_file, mdl, "INPUTDIR", inputdir, default=".") inputdir = slasher(inputdir) + ALLOC_(CS%Kh_bg_2d(isd:ied,jsd:jed)) ; CS%Kh_bg_2d(:,:) = 0.0 call MOM_read_data(trim(inputdir)//trim(filename), 'Kh', CS%Kh_bg_2d, & G%domain, timelevel=1, scale=US%m_to_L**2*US%T_to_s) call pass_var(CS%Kh_bg_2d, G%domain) @@ -2579,8 +2577,7 @@ subroutine hor_visc_init(Time, G, GV, US, param_file, diag, CS, MEKE, ADp) cmor_field_name='dispkexyfo', & cmor_long_name='Depth integrated ocean kinetic energy dissipation due to lateral friction',& cmor_standard_name='ocean_kinetic_energy_dissipation_per_unit_area_due_to_xy_friction') - if (CS%Laplacian .or. get_all) then - endif + end subroutine hor_visc_init !> Calculates factors in the anisotropic orientation tensor to be align with the grid. @@ -2618,8 +2615,9 @@ subroutine smooth_GME(CS,G,GME_flux_h,GME_flux_q) do s=1,1 ! Update halos if (present(GME_flux_h)) then + !### Work on a wider halo to eliminate this blocking send! call pass_var(GME_flux_h, G%Domain) - GME_flux_h_original = GME_flux_h + GME_flux_h_original(:,:) = GME_flux_h(:,:) ! apply smoothing on GME do j = G%jsc, G%jec do i = G%isc, G%iec @@ -2631,6 +2629,7 @@ subroutine smooth_GME(CS,G,GME_flux_h,GME_flux_q) ws = 0.125 * G%mask2dT(i,j-1) wn = 0.125 * G%mask2dT(i,j+1) wc = 1.0 - (ww+we+wn+ws) + !### Add parentheses to make this rotationally invariant. GME_flux_h(i,j) = wc * GME_flux_h_original(i,j) & + ww * GME_flux_h_original(i-1,j) & + we * GME_flux_h_original(i+1,j) & @@ -2641,8 +2640,9 @@ subroutine smooth_GME(CS,G,GME_flux_h,GME_flux_q) endif ! Update halos if (present(GME_flux_q)) then + !### Work on a wider halo to eliminate this blocking send! call pass_var(GME_flux_q, G%Domain, position=CORNER, complete=.true.) - GME_flux_q_original = GME_flux_q + GME_flux_q_original(:,:) = GME_flux_q(:,:) ! apply smoothing on GME do J = G%JscB, G%JecB do I = G%IscB, G%IecB @@ -2654,6 +2654,7 @@ subroutine smooth_GME(CS,G,GME_flux_h,GME_flux_q) ws = 0.125 * G%mask2dBu(I,J-1) wn = 0.125 * G%mask2dBu(I,J+1) wc = 1.0 - (ww+we+wn+ws) + !### Add parentheses to make this rotationally invariant. GME_flux_q(I,J) = wc * GME_flux_q_original(I,J) & + ww * GME_flux_q_original(I-1,J) & + we * GME_flux_q_original(I+1,J) & From 4c25855e016641ea07688ceed3227bf1ae2c720f Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Thu, 28 Oct 2021 21:45:03 -0400 Subject: [PATCH 2/2] Restore temporary variables Undid changes that eliminated temporary variables to facilitate performance profiling, and restored the "Knuth" convention about the placement of the line breaks relative to "+" in these expressions. Nothing changes. --- src/parameterizations/lateral/MOM_hor_visc.F90 | 18 +++++++++++------- 1 file changed, 11 insertions(+), 7 deletions(-) diff --git a/src/parameterizations/lateral/MOM_hor_visc.F90 b/src/parameterizations/lateral/MOM_hor_visc.F90 index 35a1f511b3..db2514576d 100644 --- a/src/parameterizations/lateral/MOM_hor_visc.F90 +++ b/src/parameterizations/lateral/MOM_hor_visc.F90 @@ -338,6 +338,8 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, ! viscosity. Here set equal to nondimensional Laplacian Leith constant. ! This is set equal to zero if modified Leith is not used. real :: Shear_mag_bc ! Shear_mag value in backscatter [T-1 ~> s-1] + real :: sh_xx_sq ! Square of tension (sh_xx) [T-2 ~> s-2] + real :: sh_xy_sq ! Square of shearing strain (sh_xy) [T-2 ~> s-2] real :: h2uq, h2vq ! temporary variables [H2 ~> m2 or kg2 m-4]. real :: hu, hv ! Thicknesses interpolated by arithmetic means to corner ! points; these are first interpolated to u or v velocity @@ -571,7 +573,7 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, !$OMP grad_vel_mag_h, grad_vel_mag_q, & !$OMP grad_vel_mag_bt_h, grad_vel_mag_bt_q, grad_d2vel_mag_h, & !$OMP meke_res_fn, Shear_mag, Shear_mag_bc, vert_vort_mag, h_min, hrat_min, visc_bound_rem, & - !$OMP grid_Ah, grid_Kh, d_Del2u, d_Del2v, d_str, & + !$OMP sh_xx_sq, sh_xy_sq, grid_Ah, grid_Kh, d_Del2u, d_Del2v, d_str, & !$OMP Kh, Ah, AhSm, AhLth, local_strain, Sh_F_pow, & !$OMP dDel2vdx, dDel2udy, DY_dxCv, DX_dyCu, Del2vort_q, Del2vort_h, KE, & !$OMP h2uq, h2vq, hu, hv, hq, FatH, RoScl, GME_coeff & @@ -900,9 +902,10 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, if ((CS%Smagorinsky_Kh) .or. (CS%Smagorinsky_Ah)) then do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 - Shear_mag(i,j) = sqrt( sh_xx(i,j)**2 + & - 0.25*((sh_xy(I-1,J-1)**2 + sh_xy(I,J)**2) + & - (sh_xy(I-1,J)**2 + sh_xy(I,J-1)**2)) ) + sh_xx_sq = sh_xx(i,j)**2 + sh_xy_sq = 0.25 * ( (sh_xy(I-1,J-1)**2 + sh_xy(I,J)**2) & + + (sh_xy(I-1,J)**2 + sh_xy(I,J-1)**2) ) + Shear_mag(i,j) = sqrt(sh_xx_sq + sh_xy_sq) enddo ; enddo endif @@ -1179,9 +1182,10 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, if ((CS%Smagorinsky_Kh) .or. (CS%Smagorinsky_Ah)) then do J=js-1,Jeq ; do I=is-1,Ieq - Shear_mag(I,J) = sqrt( sh_xy(I,J)**2 + & - 0.25*((sh_xx(i,j)**2 + sh_xx(i+1,j+1)**2) + & - (sh_xx(i,j+1)**2 + sh_xx(i+1,j)**2)) ) + sh_xy_sq = sh_xy(I,J)**2 + sh_xx_sq = 0.25 * ( (sh_xx(i,j)**2 + sh_xx(i+1,j+1)**2) & + + (sh_xx(i,j+1)**2 + sh_xx(i+1,j)**2) ) + Shear_mag(I,J) = sqrt(sh_xy_sq + sh_xx_sq) enddo ; enddo endif