From 32f631623eccfc88a4d3dd0c97001dc5eeef4cc2 Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Thu, 29 Feb 2024 17:09:31 -0500 Subject: [PATCH 1/7] +Rotationally symmetric neutral diffusion option Added the option to do neutral tracer diffusion with expressions that satisfy rotational symmetry. This option is enabled by setting the new runtime parameter NDIFF_ANSWER_DATE to be greater than 20240330. By default, this parameter is set to use the previous expressions, although the default may be changed later to follow DEFAULT_ANSWER_DATE. By default all answers are bitwise identical, but there are changes to some MOM_parameter_doc files due to the introduction of a new runtime parameter. --- src/tracer/MOM_neutral_diffusion.F90 | 110 ++++++++++++++++++++------- 1 file changed, 84 insertions(+), 26 deletions(-) diff --git a/src/tracer/MOM_neutral_diffusion.F90 b/src/tracer/MOM_neutral_diffusion.F90 index 87a8881b10..b64c665c87 100644 --- a/src/tracer/MOM_neutral_diffusion.F90 +++ b/src/tracer/MOM_neutral_diffusion.F90 @@ -119,6 +119,10 @@ module MOM_neutral_diffusion !! for remapping. Values below 20190101 recover the remapping !! answers from 2018, while higher values use more robust !! forms of the same remapping expressions. + integer :: ndiff_answer_date !< The vintage of the order of arithmetic to use for the neutral + !! diffusion. Values of 20240330 or below recover the answers + !! from the original form of this code, while higher values use + !! mathematically equivalent expressions that recover rotational symmetry. type(KPP_CS), pointer :: KPP_CSp => NULL() !< KPP control structure needed to get BLD type(energetic_PBL_CS), pointer :: energetic_PBL_CSp => NULL()!< ePBL control structure needed to get MLD end type neutral_diffusion_CS @@ -200,6 +204,16 @@ logical function neutral_diffusion_init(Time, G, GV, US, param_file, diag, EOS, "transports that were unmasked, as used prior to Jan 2018. This is not "//& "recommended.", default=.false.) + call get_param(param_file, mdl, "DEFAULT_ANSWER_DATE", default_answer_date, & + "This sets the default value for the various _ANSWER_DATE parameters.", & + default=99991231) + call get_param(param_file, mdl, "NDIFF_ANSWER_DATE", CS%ndiff_answer_date, & + "The vintage of the order of arithmetic to use for the neutral diffusion. "//& + "Values of 20240330 or below recover the answers from the original form of the "//& + "neutral diffusion code, while higher values use mathematically equivalent "//& + "expressions that recover rotational symmetry.", & + default=20240101) !### Change this default later to default_answer_date. + ! Initialize and configure remapping if ( .not.CS%continuous_reconstruction ) then call get_param(param_file, mdl, "NDIFF_BOUNDARY_EXTRAP", boundary_extrap, & @@ -211,9 +225,6 @@ logical function neutral_diffusion_init(Time, G, GV, US, param_file, diag, EOS, "for vertical remapping for all variables. "//& "It can be one of the following schemes: "//& trim(remappingSchemesDoc), default=remappingDefaultScheme) - call get_param(param_file, mdl, "DEFAULT_ANSWER_DATE", default_answer_date, & - "This sets the default value for the various _ANSWER_DATE parameters.", & - default=99991231) call get_param(param_file, mdl, "REMAPPING_ANSWER_DATE", CS%remap_answer_date, & "The vintage of the expressions and order of arithmetic to use for remapping. "//& "Values below 20190101 result in the use of older, less accurate expressions "//& @@ -623,6 +634,18 @@ subroutine neutral_diffusion(G, GV, h, Coef_x, Coef_y, dt, Reg, US, CS) real, dimension(SZK_(GV)) :: dTracer ! Change in tracer concentration due to neutral diffusion ! [H L2 conc ~> m3 conc or kg conc]. For temperature ! these units are [C H L2 ~> degC m3 or degC kg]. + real, dimension(SZK_(GV)) :: dTracer_N ! Change in tracer concentration due to neutral diffusion + ! into a cell via its logically northern face, in + ! [H L2 conc ~> m3 conc or kg conc]. + real, dimension(SZK_(GV)) :: dTracer_S ! Change in tracer concentration due to neutral diffusion + ! into a cell via its logically southern face, in + ! [H L2 conc ~> m3 conc or kg conc]. + real, dimension(SZK_(GV)) :: dTracer_E ! Change in tracer concentration due to neutral diffusion + ! into a cell via its logically eastern face, in + ! [H L2 conc ~> m3 conc or kg conc]. + real, dimension(SZK_(GV)) :: dTracer_W ! Change in tracer concentration due to neutral diffusion + ! into a cell via its logically western face, in + ! [H L2 conc ~> m3 conc or kg conc]. real :: normalize ! normalization used for averaging Coef_x and Coef_y to t-points [nondim]. type(tracer_type), pointer :: Tracer => NULL() ! Pointer to the current tracer @@ -800,21 +823,39 @@ subroutine neutral_diffusion(G, GV, h, Coef_x, Coef_y, dt, Reg, US, CS) endif endif - ! Update the tracer concentration from divergence of neutral diffusive flux components + ! Update the tracer concentration from divergence of neutral diffusive flux components, noting + ! that uFlx and vFlx use an unexpected sign convention. if (CS%KhTh_use_ebt_struct) then do j = G%jsc,G%jec ; do i = G%isc,G%iec if (G%mask2dT(i,j)>0.) then - dTracer(:) = 0. - do ks = 1,CS%nsurf-1 - k = CS%uKoL(I,j,ks) - dTracer(k) = dTracer(k) + uFlx(I,j,ks) - k = CS%uKoR(I-1,j,ks) - dTracer(k) = dTracer(k) - uFlx(I-1,j,ks) - k = CS%vKoL(i,J,ks) - dTracer(k) = dTracer(k) + vFlx(i,J,ks) - k = CS%vKoR(i,J-1,ks) - dTracer(k) = dTracer(k) - vFlx(i,J-1,ks) - enddo + if (CS%ndiff_answer_date <= 20240330) then + dTracer(:) = 0. + do ks = 1,CS%nsurf-1 + k = CS%uKoL(I,j,ks) + dTracer(k) = dTracer(k) + uFlx(I,j,ks) + k = CS%uKoR(I-1,j,ks) + dTracer(k) = dTracer(k) - uFlx(I-1,j,ks) + k = CS%vKoL(i,J,ks) + dTracer(k) = dTracer(k) + vFlx(i,J,ks) + k = CS%vKoR(i,J-1,ks) + dTracer(k) = dTracer(k) - vFlx(i,J-1,ks) + enddo + else ! This form recovers rotational symmetry. + dTracer_N(:) = 0.0 ; dTracer_S(:) = 0.0 ; dTracer_E(:) = 0.0 ; dTracer_W(:) = 0.0 + do ks = 1,CS%nsurf-1 + k = CS%uKoL(I,j,ks) + dTracer_E(k) = dTracer_E(k) + uFlx(I,j,ks) + k = CS%uKoR(I-1,j,ks) + dTracer_W(k) = dTracer_W(k) - uFlx(I-1,j,ks) + k = CS%vKoL(i,J,ks) + dTracer_N(k) = dTracer_N(k) + vFlx(i,J,ks) + k = CS%vKoR(i,J-1,ks) + dTracer_S(k) = dTracer_S(k) - vFlx(i,J-1,ks) + enddo + do k = 1, GV%ke + dTracer(k) = (dTracer_N(k) + dTracer_S(k)) + (dTracer_E(k) + dTracer_W(k)) + enddo + endif do k = 1, GV%ke tracer%t(i,j,k) = tracer%t(i,j,k) + dTracer(k) * & ( G%IareaT(i,j) / ( h(i,j,k) + GV%H_subroundoff ) ) @@ -832,17 +873,34 @@ subroutine neutral_diffusion(G, GV, h, Coef_x, Coef_y, dt, Reg, US, CS) else do j = G%jsc,G%jec ; do i = G%isc,G%iec if (G%mask2dT(i,j)>0.) then - dTracer(:) = 0. - do ks = 1,CS%nsurf-1 - k = CS%uKoL(I,j,ks) - dTracer(k) = dTracer(k) + Coef_x(I,j,1) * uFlx(I,j,ks) - k = CS%uKoR(I-1,j,ks) - dTracer(k) = dTracer(k) - Coef_x(I-1,j,1) * uFlx(I-1,j,ks) - k = CS%vKoL(i,J,ks) - dTracer(k) = dTracer(k) + Coef_y(i,J,1) * vFlx(i,J,ks) - k = CS%vKoR(i,J-1,ks) - dTracer(k) = dTracer(k) - Coef_y(i,J-1,1) * vFlx(i,J-1,ks) - enddo + if (CS%ndiff_answer_date <= 20240330) then + dTracer(:) = 0. + do ks = 1,CS%nsurf-1 + k = CS%uKoL(I,j,ks) + dTracer(k) = dTracer(k) + Coef_x(I,j,1) * uFlx(I,j,ks) + k = CS%uKoR(I-1,j,ks) + dTracer(k) = dTracer(k) - Coef_x(I-1,j,1) * uFlx(I-1,j,ks) + k = CS%vKoL(i,J,ks) + dTracer(k) = dTracer(k) + Coef_y(i,J,1) * vFlx(i,J,ks) + k = CS%vKoR(i,J-1,ks) + dTracer(k) = dTracer(k) - Coef_y(i,J-1,1) * vFlx(i,J-1,ks) + enddo + else ! This form recovers rotational symmetry. + dTracer_N(:) = 0.0 ; dTracer_S(:) = 0.0 ; dTracer_E(:) = 0.0 ; dTracer_W(:) = 0.0 + do ks = 1,CS%nsurf-1 + k = CS%uKoL(I,j,ks) + dTracer_E(k) = dTracer_E(k) + Coef_x(I,j,1) * uFlx(I,j,ks) + k = CS%uKoR(I-1,j,ks) + dTracer_W(k) = dTracer_W(k) - Coef_x(I-1,j,1) * uFlx(I-1,j,ks) + k = CS%vKoL(i,J,ks) + dTracer_N(k) = dTracer_N(k) + Coef_y(i,J,1) * vFlx(i,J,ks) + k = CS%vKoR(i,J-1,ks) + dTracer_S(k) = dTracer_S(k) - Coef_y(i,J-1,1) * vFlx(i,J-1,ks) + enddo + do k = 1, GV%ke + dTracer(k) = (dTracer_N(k) + dTracer_S(k)) + (dTracer_E(k) + dTracer_W(k)) + enddo + endif do k = 1, GV%ke tracer%t(i,j,k) = tracer%t(i,j,k) + dTracer(k) * & ( G%IareaT(i,j) / ( h(i,j,k) + GV%H_subroundoff ) ) From 812510bc4906426ed2f0a2d11cac1629e24369f0 Mon Sep 17 00:00:00 2001 From: He Wang Date: Mon, 4 Mar 2024 00:01:57 -0500 Subject: [PATCH 2/7] Add an option in hor_visc to use cont thickness Runtime parameter USE_CONT_THICKNESS is added in hor_visc to let it use velocity-point thickness consistent with the continuity solver. Thicknesses are borrowed from BT_cont so only split mode is supported. --- src/core/MOM_dynamics_split_RK2.F90 | 5 ++-- .../lateral/MOM_hor_visc.F90 | 24 +++++++++++++++++-- 2 files changed, 25 insertions(+), 4 deletions(-) diff --git a/src/core/MOM_dynamics_split_RK2.F90 b/src/core/MOM_dynamics_split_RK2.F90 index 11b1eb5c16..0557ec7cd5 100644 --- a/src/core/MOM_dynamics_split_RK2.F90 +++ b/src/core/MOM_dynamics_split_RK2.F90 @@ -851,7 +851,7 @@ subroutine step_MOM_dyn_split_RK2(u_inst, v_inst, h, tv, visc, Time_local, dt, f call horizontal_viscosity(u_av, v_av, h_av, CS%diffu, CS%diffv, & MEKE, Varmix, G, GV, US, CS%hor_visc, & OBC=CS%OBC, BT=CS%barotropic_CSp, TD=thickness_diffuse_CSp, & - ADp=CS%ADp) + ADp=CS%ADp, hu_cont=CS%BT_cont%h_u, hv_cont=CS%BT_cont%h_v) call cpu_clock_end(id_clock_horvisc) if (showCallTree) call callTree_wayPoint("done with horizontal_viscosity (step_MOM_dyn_split_RK2)") @@ -1518,7 +1518,8 @@ subroutine initialize_dyn_split_RK2(u, v, h, uh, vh, eta, Time, G, GV, US, param if (.not. query_initialized(CS%diffu, "diffu", restart_CS) .or. & .not. query_initialized(CS%diffv, "diffv", restart_CS)) then call horizontal_viscosity(u, v, h, CS%diffu, CS%diffv, MEKE, VarMix, G, GV, US, CS%hor_visc, & - OBC=CS%OBC, BT=CS%barotropic_CSp, TD=thickness_diffuse_CSp) + OBC=CS%OBC, BT=CS%barotropic_CSp, TD=thickness_diffuse_CSp, & + hu_cont=CS%BT_cont%h_u, hv_cont=CS%BT_cont%h_v) call set_initialized(CS%diffu, "diffu", restart_CS) call set_initialized(CS%diffv, "diffv", restart_CS) endif diff --git a/src/parameterizations/lateral/MOM_hor_visc.F90 b/src/parameterizations/lateral/MOM_hor_visc.F90 index e3249afb73..96b8bd2a9e 100644 --- a/src/parameterizations/lateral/MOM_hor_visc.F90 +++ b/src/parameterizations/lateral/MOM_hor_visc.F90 @@ -111,7 +111,7 @@ module MOM_hor_visc !! 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] - + logical :: use_cont_thick !< If true, thickness at velocity points adopts h[uv] in BT_cont from continuity solver. type(ZB2020_CS) :: ZB2020 !< Zanna-Bolton 2020 control structure. logical :: use_ZB2020 !< If true, use Zanna-Bolton 2020 parameterization. @@ -239,7 +239,7 @@ module MOM_hor_visc !! v[is-2:ie+2,js-2:je+2] !! h[is-1:ie+1,js-1:je+1] subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, & - CS, OBC, BT, TD, ADp) + CS, OBC, BT, TD, ADp, hu_cont, hv_cont) type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure. type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure. real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), & @@ -263,6 +263,10 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, type(barotropic_CS), intent(in), optional :: BT !< Barotropic control structure type(thickness_diffuse_CS), intent(in), optional :: TD !< Thickness diffusion control structure type(accel_diag_ptrs), intent(in), optional :: ADp !< Acceleration diagnostics + real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), & + optional, intent(in) :: hu_cont !< Layer thickness at u-points [H ~> m or kg m-2]. + real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), & + optional, intent(in) :: hv_cont !< Layer thickness at v-points [H ~> m or kg m-2]. ! Local variables real, dimension(SZIB_(G),SZJ_(G)) :: & @@ -391,6 +395,7 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, logical :: apply_OBC = .false. logical :: use_MEKE_Ku logical :: use_MEKE_Au + logical :: use_cont_huv integer :: is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz integer :: i, j, k, n real :: inv_PI3, inv_PI2, inv_PI6 ! Powers of the inverse of pi [nondim] @@ -445,6 +450,8 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, use_MEKE_Ku = allocated(MEKE%Ku) use_MEKE_Au = allocated(MEKE%Au) + use_cont_huv = CS%use_cont_thick .and. present(hu_cont) .and. present(hv_cont) + rescale_Kh = .false. if (VarMix%use_variable_mixing) then rescale_Kh = VarMix%Resoln_scaled_Kh @@ -658,6 +665,16 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, enddo ; enddo endif + ! The following should obviously be combined with the previous block if adopted. + if (use_cont_huv) then + do j=js-2,je+2 ; do I=Isq-1,Ieq+1 + h_u(I,j) = hu_cont(I,j,k) + enddo ; enddo + do J=Jsq-1,Jeq+1 ; do i=is-2,ie+2 + h_v(i,J) = hv_cont(i,J,k) + enddo ; enddo + endif + ! Adjust contributions to shearing strain and interpolated values of ! thicknesses on open boundaries. if (apply_OBC) then ; do n=1,OBC%number_of_segments @@ -1969,6 +1986,9 @@ subroutine hor_visc_init(Time, G, GV, US, param_file, diag, CS, ADp) if (.not.GV%Boussinesq) CS%answer_date = max(CS%answer_date, 20230701) call get_param(param_file, mdl, "DEBUG", CS%debug, default=.false.) + call get_param(param_file, mdl, "USE_CONT_THICKNESS", CS%use_cont_thick, & + "If true, use thickness at velocity points from continuity solver. This option"//& + "currently only works with split mode.", default=.false.) call get_param(param_file, mdl, "LAPLACIAN", CS%Laplacian, & "If true, use a Laplacian horizontal viscosity.", & default=.false.) From 37ff301a65817e7212cd7743c66fb406fc2b9223 Mon Sep 17 00:00:00 2001 From: He Wang Date: Mon, 4 Mar 2024 08:44:52 -0500 Subject: [PATCH 3/7] fix opemmp parallel --- src/parameterizations/lateral/MOM_hor_visc.F90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/parameterizations/lateral/MOM_hor_visc.F90 b/src/parameterizations/lateral/MOM_hor_visc.F90 index 96b8bd2a9e..6f707f9e87 100644 --- a/src/parameterizations/lateral/MOM_hor_visc.F90 +++ b/src/parameterizations/lateral/MOM_hor_visc.F90 @@ -561,12 +561,12 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, !$OMP CS, G, GV, US, OBC, VarMix, MEKE, u, v, h, & !$OMP is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz, & !$OMP apply_OBC, rescale_Kh, legacy_bound, find_FrictWork, & - !$OMP use_MEKE_Ku, use_MEKE_Au, & + !$OMP use_MEKE_Ku, use_MEKE_Au, use_cont_huv, & !$OMP backscat_subround, GME_effic_h, GME_effic_q, & !$OMP h_neglect, h_neglect3, inv_PI3, inv_PI6, & !$OMP diffu, diffv, Kh_h, Kh_q, Ah_h, Ah_q, FrictWork, FrictWork_GME, & !$OMP div_xx_h, sh_xx_h, vort_xy_q, sh_xy_q, GME_coeff_h, GME_coeff_q, & - !$OMP KH_u_GME, KH_v_GME, grid_Re_Kh, grid_Re_Ah, NoSt, ShSt & + !$OMP KH_u_GME, KH_v_GME, grid_Re_Kh, grid_Re_Ah, NoSt, ShSt, hu_cont, hv_cont & !$OMP ) & !$OMP private( & !$OMP i, j, k, n, & From 2b59089ea5561b14cb55c092ec78b1c0e1864a2c Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Tue, 5 Mar 2024 05:52:47 -0500 Subject: [PATCH 4/7] Set up auxiliary domain for unrotated grid When ROTATE_INDEX is true, the model keeps both its rotated internal grid and an unrotated grid for setting initialization and forcing. When the model is run in symmetric memory mode, there is an auxiliary non-symmetric domain that is needed for reading in fields at velocity points. The auxiliary domain in the unrotated grid (G_in) was not being set previously, causing the model to fail to run some cases when ROTATE_INDEX was true. That auxiliary domain in the unrotated grid is now being set up properly. All answers and output are bitwise identical for any cases that worked previously. --- src/core/MOM.F90 | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/src/core/MOM.F90 b/src/core/MOM.F90 index 094e534312..595b730e0f 100644 --- a/src/core/MOM.F90 +++ b/src/core/MOM.F90 @@ -2904,6 +2904,11 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, & if (CS%rotate_index) then G_in%ke = GV%ke + ! Allocate the auxiliary non-symmetric domain for debugging or I/O purposes. + if (CS%debug .or. G_in%symmetric) then + call clone_MOM_domain(G_in%Domain, G_in%Domain_aux, symmetric=.false.) + else ; G_in%Domain_aux => G_in%Domain ; endif + allocate(u_in(G_in%IsdB:G_in%IedB, G_in%jsd:G_in%jed, nz), source=0.0) allocate(v_in(G_in%isd:G_in%ied, G_in%JsdB:G_in%JedB, nz), source=0.0) allocate(h_in(G_in%isd:G_in%ied, G_in%jsd:G_in%jed, nz), source=GV%Angstrom_H) From 5e34f486b0b50d4f7bdbf752456b355824d6e885 Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Wed, 6 Mar 2024 15:31:51 -0500 Subject: [PATCH 5/7] Set flags properly for rotated tripolar grids Added code to tht two versions of clone_MD_to_MD in the FMS1 and FMS2 versions of MOM_domain_infra.F90 to properly change the flags when a tripolar grid is rotated, so that it does not lead to misleadingly incorrect answers. With the current versions of FMS, doing grid rotation testing with a tripolar grid will lead to error messages about the incomplete implementation of FMS, but these failures are preferable to the model silently working incorrectly. All answers are bitwise identical in cases that worked correctly before. --- config_src/infra/FMS1/MOM_domain_infra.F90 | 26 ++++++++++++++++++++-- config_src/infra/FMS2/MOM_domain_infra.F90 | 26 ++++++++++++++++++++-- 2 files changed, 48 insertions(+), 4 deletions(-) diff --git a/config_src/infra/FMS1/MOM_domain_infra.F90 b/config_src/infra/FMS1/MOM_domain_infra.F90 index 2c97a0bb31..8040be27d9 100644 --- a/config_src/infra/FMS1/MOM_domain_infra.F90 +++ b/config_src/infra/FMS1/MOM_domain_infra.F90 @@ -19,7 +19,8 @@ module MOM_domain_infra use mpp_domains_mod, only : mpp_compute_block_extent use mpp_domains_mod, only : mpp_broadcast_domain, mpp_redistribute, mpp_global_field use mpp_domains_mod, only : AGRID, BGRID_NE, CGRID_NE, SCALAR_PAIR, BITWISE_EXACT_SUM -use mpp_domains_mod, only : CYCLIC_GLOBAL_DOMAIN, FOLD_NORTH_EDGE +use mpp_domains_mod, only : CYCLIC_GLOBAL_DOMAIN +use mpp_domains_mod, only : FOLD_NORTH_EDGE, FOLD_SOUTH_EDGE, FOLD_EAST_EDGE, FOLD_WEST_EDGE use mpp_domains_mod, only : To_East => WUPDATE, To_West => EUPDATE, Omit_Corners => EDGEUPDATE use mpp_domains_mod, only : To_North => SUPDATE, To_South => NUPDATE use mpp_domains_mod, only : CENTER, CORNER, NORTH_FACE => NORTH, EAST_FACE => EAST @@ -1553,6 +1554,19 @@ subroutine clone_MD_to_MD(MD_in, MOM_dom, min_halo, halo_size, symmetric, domain call get_layout_extents(MD_in, exnj, exni) MOM_dom%X_FLAGS = MD_in%Y_FLAGS ; MOM_dom%Y_FLAGS = MD_in%X_FLAGS + ! Correct the position of a tripolar grid, assuming that flags are not additive. + if (qturns == 1) then + if (MD_in%Y_FLAGS == FOLD_NORTH_EDGE) MOM_dom%X_FLAGS = FOLD_EAST_EDGE + if (MD_in%Y_FLAGS == FOLD_SOUTH_EDGE) MOM_dom%X_FLAGS = FOLD_WEST_EDGE + if (MD_in%X_FLAGS == FOLD_EAST_EDGE) MOM_dom%Y_FLAGS = FOLD_SOUTH_EDGE + if (MD_in%X_FLAGS == FOLD_WEST_EDGE) MOM_dom%Y_FLAGS = FOLD_NORTH_EDGE + elseif (qturns == 3) then + if (MD_in%Y_FLAGS == FOLD_NORTH_EDGE) MOM_dom%X_FLAGS = FOLD_WEST_EDGE + if (MD_in%Y_FLAGS == FOLD_SOUTH_EDGE) MOM_dom%X_FLAGS = FOLD_EAST_EDGE + if (MD_in%X_FLAGS == FOLD_EAST_EDGE) MOM_dom%Y_FLAGS = FOLD_NORTH_EDGE + if (MD_in%X_FLAGS == FOLD_WEST_EDGE) MOM_dom%Y_FLAGS = FOLD_SOUTH_EDGE + endif + MOM_dom%layout(:) = MD_in%layout(2:1:-1) MOM_dom%io_layout(:) = io_layout_in(2:1:-1) else @@ -1561,11 +1575,19 @@ subroutine clone_MD_to_MD(MD_in, MOM_dom, min_halo, halo_size, symmetric, domain call get_layout_extents(MD_in, exni, exnj) MOM_dom%X_FLAGS = MD_in%X_FLAGS ; MOM_dom%Y_FLAGS = MD_in%Y_FLAGS + ! Correct the position of a tripolar grid, assuming that flags are not additive. + if (qturns == 2) then + if (MD_in%Y_FLAGS == FOLD_NORTH_EDGE) MOM_dom%Y_FLAGS = FOLD_SOUTH_EDGE + if (MD_in%Y_FLAGS == FOLD_SOUTH_EDGE) MOM_dom%Y_FLAGS = FOLD_NORTH_EDGE + if (MD_in%X_FLAGS == FOLD_EAST_EDGE) MOM_dom%X_FLAGS = FOLD_WEST_EDGE + if (MD_in%X_FLAGS == FOLD_WEST_EDGE) MOM_dom%X_FLAGS = FOLD_EAST_EDGE + endif + MOM_dom%layout(:) = MD_in%layout(:) MOM_dom%io_layout(:) = io_layout_in(:) endif - ! Ensure that the points per processor are the same on the source and densitation grids. + ! Ensure that the points per processor are the same on the source and destination grids. select case (qturns) case (1) ; call invert(exni) case (2) ; call invert(exni) ; call invert(exnj) diff --git a/config_src/infra/FMS2/MOM_domain_infra.F90 b/config_src/infra/FMS2/MOM_domain_infra.F90 index ff1d888c47..76d9469e3c 100644 --- a/config_src/infra/FMS2/MOM_domain_infra.F90 +++ b/config_src/infra/FMS2/MOM_domain_infra.F90 @@ -19,7 +19,8 @@ module MOM_domain_infra use mpp_domains_mod, only : mpp_compute_block_extent use mpp_domains_mod, only : mpp_broadcast_domain, mpp_redistribute, mpp_global_field use mpp_domains_mod, only : AGRID, BGRID_NE, CGRID_NE, SCALAR_PAIR, BITWISE_EXACT_SUM -use mpp_domains_mod, only : CYCLIC_GLOBAL_DOMAIN, FOLD_NORTH_EDGE +use mpp_domains_mod, only : CYCLIC_GLOBAL_DOMAIN +use mpp_domains_mod, only : FOLD_NORTH_EDGE, FOLD_SOUTH_EDGE, FOLD_EAST_EDGE, FOLD_WEST_EDGE use mpp_domains_mod, only : To_East => WUPDATE, To_West => EUPDATE, Omit_Corners => EDGEUPDATE use mpp_domains_mod, only : To_North => SUPDATE, To_South => NUPDATE use mpp_domains_mod, only : CENTER, CORNER, NORTH_FACE => NORTH, EAST_FACE => EAST @@ -1555,6 +1556,19 @@ subroutine clone_MD_to_MD(MD_in, MOM_dom, min_halo, halo_size, symmetric, domain call get_layout_extents(MD_in, exnj, exni) MOM_dom%X_FLAGS = MD_in%Y_FLAGS ; MOM_dom%Y_FLAGS = MD_in%X_FLAGS + ! Correct the position of a tripolar grid, assuming that flags are not additive. + if (modulo(qturns, 4) == 1) then + if (MD_in%Y_FLAGS == FOLD_NORTH_EDGE) MOM_dom%X_FLAGS = FOLD_EAST_EDGE + if (MD_in%Y_FLAGS == FOLD_SOUTH_EDGE) MOM_dom%X_FLAGS = FOLD_WEST_EDGE + if (MD_in%X_FLAGS == FOLD_EAST_EDGE) MOM_dom%Y_FLAGS = FOLD_SOUTH_EDGE + if (MD_in%X_FLAGS == FOLD_WEST_EDGE) MOM_dom%Y_FLAGS = FOLD_NORTH_EDGE + elseif (modulo(qturns, 4) == 3) then + if (MD_in%Y_FLAGS == FOLD_NORTH_EDGE) MOM_dom%X_FLAGS = FOLD_WEST_EDGE + if (MD_in%Y_FLAGS == FOLD_SOUTH_EDGE) MOM_dom%X_FLAGS = FOLD_EAST_EDGE + if (MD_in%X_FLAGS == FOLD_EAST_EDGE) MOM_dom%Y_FLAGS = FOLD_NORTH_EDGE + if (MD_in%X_FLAGS == FOLD_WEST_EDGE) MOM_dom%Y_FLAGS = FOLD_SOUTH_EDGE + endif + MOM_dom%layout(:) = MD_in%layout(2:1:-1) MOM_dom%io_layout(:) = io_layout_in(2:1:-1) else @@ -1563,11 +1577,19 @@ subroutine clone_MD_to_MD(MD_in, MOM_dom, min_halo, halo_size, symmetric, domain call get_layout_extents(MD_in, exni, exnj) MOM_dom%X_FLAGS = MD_in%X_FLAGS ; MOM_dom%Y_FLAGS = MD_in%Y_FLAGS + ! Correct the position of a tripolar grid, assuming that flags are not additive. + if (modulo(qturns, 4) == 2) then + if (MD_in%Y_FLAGS == FOLD_NORTH_EDGE) MOM_dom%Y_FLAGS = FOLD_SOUTH_EDGE + if (MD_in%Y_FLAGS == FOLD_SOUTH_EDGE) MOM_dom%Y_FLAGS = FOLD_NORTH_EDGE + if (MD_in%X_FLAGS == FOLD_EAST_EDGE) MOM_dom%X_FLAGS = FOLD_WEST_EDGE + if (MD_in%X_FLAGS == FOLD_WEST_EDGE) MOM_dom%X_FLAGS = FOLD_EAST_EDGE + endif + MOM_dom%layout(:) = MD_in%layout(:) MOM_dom%io_layout(:) = io_layout_in(:) endif - ! Ensure that the points per processor are the same on the source and densitation grids. + ! Ensure that the points per processor are the same on the source and destination grids. select case (qturns) case (1) ; call invert(exni) case (2) ; call invert(exni) ; call invert(exnj) From f9372f3d66392df199ae7541c62a99f467971285 Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Sun, 10 Mar 2024 22:46:41 -0400 Subject: [PATCH 6/7] +Add get_netcdf_filename for a get_field_nc error Add get_netcdf_filename and use it to add useful details (the field and filenames in question) to a fatal error message in get_field_nc. All answers are bitwise identical, but there is a new public interface and some output is changed in cases where get_field_nc is failing. --- src/framework/MOM_io_file.F90 | 6 ++++-- src/framework/MOM_netcdf.F90 | 9 +++++++++ 2 files changed, 13 insertions(+), 2 deletions(-) diff --git a/src/framework/MOM_io_file.F90 b/src/framework/MOM_io_file.F90 index c6d86a008b..6697e56f68 100644 --- a/src/framework/MOM_io_file.F90 +++ b/src/framework/MOM_io_file.F90 @@ -34,6 +34,7 @@ module MOM_io_file use MOM_netcdf, only : write_netcdf_attribute use MOM_netcdf, only : get_netcdf_size use MOM_netcdf, only : get_netcdf_fields +use MOM_netcdf, only : get_netcdf_filename use MOM_netcdf, only : read_netcdf_field use MOM_error_handler, only : MOM_error, FATAL @@ -1757,8 +1758,9 @@ subroutine get_field_nc(handle, label, values, rescale) ! NOTE: Data on face and vertex points is not yet supported. This is a ! temporary check to detect such cases, but may be removed in the future. if (.not. (compute_domain .or. data_domain)) & - call MOM_error(FATAL, 'get_field_nc: Only compute and data domains ' // & - 'are currently supported.') + call MOM_error(FATAL, 'get_field_nc trying to read '//trim(label)//' from '//& + trim(get_netcdf_filename(handle%handle_nc))//& + ': Only compute and data domains are currently supported.') field_nc = handle%fields%get(label) diff --git a/src/framework/MOM_netcdf.F90 b/src/framework/MOM_netcdf.F90 index 95e6aa7bb7..8d6534e5dd 100644 --- a/src/framework/MOM_netcdf.F90 +++ b/src/framework/MOM_netcdf.F90 @@ -39,6 +39,7 @@ module MOM_netcdf public :: write_netcdf_attribute public :: get_netcdf_size public :: get_netcdf_fields +public :: get_netcdf_filename public :: read_netcdf_field @@ -722,6 +723,14 @@ subroutine get_netcdf_fields(handle, axes, fields) fields(:) = vars(:nfields) end subroutine get_netcdf_fields +!> Return the name of a file from a netCDF handle +function get_netcdf_filename(handle) + type(netcdf_file_type), intent(in) :: handle !< A netCDF file handle + character(len=:), allocatable :: get_netcdf_filename !< The name of the file that this handle refers to. + + get_netcdf_filename = handle%filename + +end function !> Read the values of a field from a netCDF file subroutine read_netcdf_field(handle, field, values, bounds) From a3fd1f3a7466e886b45af116ae93db6111e4c644 Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Mon, 18 Mar 2024 13:15:08 -0400 Subject: [PATCH 7/7] +Rotationally symmetric epipycnal diffusion option Added the option to do epipycnal tracer diffusion between a bulk mixed layer and the interior ocean with expressions that satisfy rotational symmetry. This option is enabled by setting the new runtime parameter HOR_DIFF_ANSWER_DATE to be greater than 20240330. By default, this parameter is set to use the previous expressions, although the default may be changed later to follow DEFAULT_ANSWER_DATE. Also corrected two bugs with the tracer limits used to repartition the fluxes between layers within tracer_epipycnal_ML_diff; this correction is enables by setting the new HOR_DIFF_LIMIT_BUG to .false., but to retain previous answers by default it is set to true. By default all answers are bitwise identical, but there are changes to some MOM_parameter_doc files due to the introduction of two new runtime parameters. --- src/tracer/MOM_tracer_hor_diff.F90 | 217 ++++++++++++++++++++--------- 1 file changed, 149 insertions(+), 68 deletions(-) diff --git a/src/tracer/MOM_tracer_hor_diff.F90 b/src/tracer/MOM_tracer_hor_diff.F90 index 6f4e5d0f90..732a42e44b 100644 --- a/src/tracer/MOM_tracer_hor_diff.F90 +++ b/src/tracer/MOM_tracer_hor_diff.F90 @@ -65,6 +65,14 @@ module MOM_tracer_hor_diff !! tracer_hor_diff. logical :: recalc_neutral_surf !< If true, recalculate the neutral surfaces if CFL has been !! exceeded + logical :: limit_bug !< If true and the answer date is 20240330 or below, use a + !! rotational symmetry breaking bug when limiting the tracer + !! properties in tracer_epipycnal_ML_diff. + integer :: answer_date !< The vintage of the order of arithmetic to use for the tracer + !! diffusion. Values of 20240330 or below recover the answers + !! from the original form of this code, while higher values use + !! mathematically equivalent expressions that recover rotational symmetry + !! when DIFFUSE_ML_TO_INTERIOR is true. type(neutral_diffusion_CS), pointer :: neutral_diffusion_CSp => NULL() !< Control structure for neutral diffusion. type(hbd_CS), pointer :: hor_bnd_diffusion_CSp => NULL() !< Control structure for !! horizontal boundary diffusion. @@ -678,7 +686,7 @@ subroutine tracer_epipycnal_ML_diff(h, dt, Tr, ntr, khdt_epi_x, khdt_epi_y, G, & real, dimension(SZI_(G),SZJB_(G)), intent(in) :: khdt_epi_y !< Meridional epipycnal diffusivity times !! a time step and the ratio of the open face width over !! the distance between adjacent tracer points [L2 ~> m2] - type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type + type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type type(tracer_hor_diff_CS), intent(inout) :: CS !< module control structure type(thermo_var_ptrs), intent(in) :: tv !< thermodynamic structure integer, intent(in) :: num_itts !< number of iterations (usually=1) @@ -706,13 +714,16 @@ subroutine tracer_epipycnal_ML_diff(h, dt, Tr, ntr, khdt_epi_x, khdt_epi_y, G, & k0b_Lv, k0a_Lv, & ! The original k-indices of the layers that participate k0b_Rv, k0a_Rv ! in each pair of mixing at v-faces. - !### Accumulating the converge into this array one face at a time may lead to a lack of rotational symmetry. - real, dimension(SZI_(G), SZJ_(G), SZK_(GV)) :: & - tr_flux_conv ! The flux convergence of tracers [conc H L2 ~> conc m3 or conc kg] + real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: & + tr_flux_N, & ! The tracer flux through the northern face [conc H L2 ~> conc m3 or conc kg] + tr_flux_S, & ! The tracer flux through the southern face [conc H L2 ~> conc m3 or conc kg] + tr_flux_E, & ! The tracer flux through the eastern face [conc H L2 ~> conc m3 or conc kg] + tr_flux_W, & ! The tracer flux through the western face [conc H L2 ~> conc m3 or conc kg] + tr_flux_conv ! The flux convergence of tracers [conc H L2 ~> conc m3 or conc kg] ! The following 3-d arrays were created in 2014 in MOM6 PR#12 to facilitate openMP threading - ! on an i-loop, which might have been ill advised. The k-size extents here might also be problematic. - real, dimension(SZI_(G),SZJB_(G),SZK_(GV)) :: & + ! on an i-loop, which might have been ill advised. + real, dimension(SZI_(G),SZJB_(G),SZK_(GV)*2) :: & Tr_flux_3d, & ! The tracer flux through pairings at meridional faces [conc H L2 ~> conc m3 or conc kg] Tr_adj_vert_L, & ! Vertical adjustments to which layer the fluxes go into in the southern ! columns at meridional face [conc H L2 ~> conc m3 or conc kg] @@ -815,6 +826,7 @@ subroutine tracer_epipycnal_ML_diff(h, dt, Tr, ntr, khdt_epi_x, khdt_epi_y, G, & do k=2,nkmb ; do j=js-2,je+2 ; do i=is-2,ie+2 if (Rml_max(i,j) < rho_coord(i,j,k)) Rml_max(i,j) = rho_coord(i,j,k) enddo ; enddo ; enddo + ! Use bracketing and bisection to find the k-level that the densest of the ! mixed and buffer layer corresponds to, such that: ! GV%Rlay(max_kRho-1) < Rml_max <= GV%Rlay(max_kRho) @@ -1191,12 +1203,7 @@ subroutine tracer_epipycnal_ML_diff(h, dt, Tr, ntr, khdt_epi_x, khdt_epi_y, G, & endif ; enddo ; enddo ! i- & j- loops over meridional faces. -! The tracer-specific calculations start here. - - ! Zero out tracer tendencies. - do k=1,PEmax_kRho ; do j=js-1,je+1 ; do i=is-1,ie+1 - tr_flux_conv(i,j,k) = 0.0 - enddo ; enddo ; enddo + ! The tracer-specific calculations start here. do itt=1,max_itt @@ -1205,12 +1212,19 @@ subroutine tracer_epipycnal_ML_diff(h, dt, Tr, ntr, khdt_epi_x, khdt_epi_y, G, & endif do m=1,ntr -!$OMP parallel do default(none) shared(is,ie,js,je,G,Tr,nkmb,nPu,m,max_kRho,nz,h,h_exclude, & -!$OMP k0b_Lu,k0b_Ru,deep_wt_Lu,k0a_Lu,deep_wt_Ru,k0a_Ru, & -!$OMP hP_Lu,hP_Ru,I_maxitt,khdt_epi_x,tr_flux_conv,Idt) & -!$OMP private(Tr_min_face,Tr_max_face,kLa,kLb,kRa,kRb,Tr_La, & -!$OMP Tr_Lb,Tr_Ra,Tr_Rb,Tr_av_L,wt_b,Tr_av_R,h_L,h_R, & -!$OMP Tr_flux,Tr_adj_vert,wt_a,vol) + ! Zero out tracer tendencies. + if (CS%answer_date <= 20240330) then + tr_flux_conv(:,:,:) = 0.0 + else + tr_flux_N(:,:,:) = 0.0 ; tr_flux_S(:,:,:) = 0.0 + tr_flux_E(:,:,:) = 0.0 ; tr_flux_W(:,:,:) = 0.0 + endif + tr_flux_3d(:,:,:) = 0.0 + tr_adj_vert_R(:,:,:) = 0.0 ; tr_adj_vert_L(:,:,:) = 0.0 + + !$OMP parallel do default(shared) private(Tr_min_face,Tr_max_face,kLa,kLb,kRa,kRb,Tr_La, & + !$OMP Tr_Lb,Tr_Ra,Tr_Rb,Tr_av_L,wt_b,Tr_av_R,h_L,h_R, & + !$OMP Tr_flux,Tr_adj_vert,wt_a,vol) do j=js,je ; do I=is-1,ie ; if (G%mask2dCu(I,j) > 0.0) then ! Determine the fluxes through the zonal faces. @@ -1230,7 +1244,11 @@ subroutine tracer_epipycnal_ML_diff(h, dt, Tr, ntr, khdt_epi_x, khdt_epi_y, G, & kRb = kRa ; if (max_kRho(i+1,j) < nz) kRb = max_kRho(i+1,j)+1 Tr_La = Tr_min_face ; Tr_Lb = Tr_La ; Tr_Ra = Tr_La ; Tr_Rb = Tr_La if (h(i,j,kLa) > h_exclude) Tr_La = Tr(m)%t(i,j,kLa) - if (h(i,j,kLb) > h_exclude) Tr_La = Tr(m)%t(i,j,kLb) + if ((CS%answer_date <= 20240330) .and. CS%limit_bug) then + if (h(i,j,kLb) > h_exclude) Tr_La = Tr(m)%t(i,j,kLb) + else + if (h(i,j,kLb) > h_exclude) Tr_Lb = Tr(m)%t(i,j,kLb) + endif if (h(i+1,j,kRa) > h_exclude) Tr_Ra = Tr(m)%t(i+1,j,kRa) if (h(i+1,j,kRb) > h_exclude) Tr_Rb = Tr(m)%t(i+1,j,kRb) Tr_min_face = min(Tr_min_face, Tr_La, Tr_Lb, Tr_Ra, Tr_Rb) @@ -1264,12 +1282,20 @@ subroutine tracer_epipycnal_ML_diff(h, dt, Tr, ntr, khdt_epi_x, khdt_epi_y, G, & endif h_L = hP_Lu(j)%p(I,k) ; h_R = hP_Ru(j)%p(I,k) - Tr_flux = I_maxitt * khdt_epi_x(I,j) * (Tr_av_L - Tr_av_R) * & - ((2.0 * h_L * h_R) / (h_L + h_R)) - + if (CS%answer_date <= 20240330) then + Tr_flux = I_maxitt * khdt_epi_x(I,j) * (Tr_av_L - Tr_av_R) * & + ((2.0 * h_L * h_R) / (h_L + h_R)) + else + Tr_flux = I_maxitt * ((2.0 * h_L * h_R) / (h_L + h_R)) * & + khdt_epi_x(I,j) * (Tr_av_L - Tr_av_R) + endif if (deep_wt_Lu(j)%p(I,k) >= 1.0) then - tr_flux_conv(i,j,kLb) = tr_flux_conv(i,j,kLb) - Tr_flux + if (CS%answer_date <= 20240330) then + tr_flux_conv(i,j,kLb) = tr_flux_conv(i,j,kLb) - Tr_flux + else + tr_flux_E(i,j,kLb) = tr_flux_E(i,j,kLb) + Tr_flux + endif else Tr_adj_vert = 0.0 wt_b = deep_wt_Lu(j)%p(I,k) ; wt_a = 1.0 - wt_b @@ -1299,12 +1325,21 @@ subroutine tracer_epipycnal_ML_diff(h, dt, Tr, ntr, khdt_epi_x, khdt_epi_y, G, & endif endif - tr_flux_conv(i,j,kLa) = tr_flux_conv(i,j,kLa) - (wt_a*Tr_flux + Tr_adj_vert) - tr_flux_conv(i,j,kLb) = tr_flux_conv(i,j,kLb) - (wt_b*Tr_flux - Tr_adj_vert) + if (CS%answer_date <= 20240330) then + tr_flux_conv(i,j,kLa) = tr_flux_conv(i,j,kLa) - (wt_a*Tr_flux + Tr_adj_vert) + tr_flux_conv(i,j,kLb) = tr_flux_conv(i,j,kLb) - (wt_b*Tr_flux - Tr_adj_vert) + else + tr_flux_E(i,j,kLa) = tr_flux_E(i,j,kLa) + (wt_a*Tr_flux + Tr_adj_vert) + tr_flux_E(i,j,kLb) = tr_flux_E(i,j,kLb) + (wt_b*Tr_flux - Tr_adj_vert) + endif endif if (deep_wt_Ru(j)%p(I,k) >= 1.0) then - tr_flux_conv(i+1,j,kRb) = tr_flux_conv(i+1,j,kRb) + Tr_flux + if (CS%answer_date <= 20240330) then + tr_flux_conv(i+1,j,kRb) = tr_flux_conv(i+1,j,kRb) + Tr_flux + else + tr_flux_W(i+1,j,kRb) = tr_flux_W(i+1,j,kRb) + Tr_flux + endif else Tr_adj_vert = 0.0 wt_b = deep_wt_Ru(j)%p(I,k) ; wt_a = 1.0 - wt_b @@ -1334,23 +1369,22 @@ subroutine tracer_epipycnal_ML_diff(h, dt, Tr, ntr, khdt_epi_x, khdt_epi_y, G, & endif endif - tr_flux_conv(i+1,j,kRa) = tr_flux_conv(i+1,j,kRa) + & - (wt_a*Tr_flux - Tr_adj_vert) - tr_flux_conv(i+1,j,kRb) = tr_flux_conv(i+1,j,kRb) + & - (wt_b*Tr_flux + Tr_adj_vert) + if (CS%answer_date <= 20240330) then + tr_flux_conv(i+1,j,kRa) = tr_flux_conv(i+1,j,kRa) + (wt_a*Tr_flux - Tr_adj_vert) + tr_flux_conv(i+1,j,kRb) = tr_flux_conv(i+1,j,kRb) + (wt_b*Tr_flux + Tr_adj_vert) + else + tr_flux_W(i+1,j,kRa) = tr_flux_W(i+1,j,kRa) + (wt_a*Tr_flux - Tr_adj_vert) + tr_flux_W(i+1,j,kRb) = tr_flux_W(i+1,j,kRb) + (wt_b*Tr_flux + Tr_adj_vert) + endif endif if (associated(Tr(m)%df2d_x)) & Tr(m)%df2d_x(I,j) = Tr(m)%df2d_x(I,j) + Tr_flux * Idt enddo ! Loop over pairings at faces. endif ; enddo ; enddo ! i- & j- loops over zonal faces. -!$OMP parallel do default(none) shared(is,ie,js,je,G,Tr,nkmb,nPv,m,max_kRho,nz,h,h_exclude, & -!$OMP k0b_Lv,k0b_Rv,deep_wt_Lv,k0a_Lv,deep_wt_Rv,k0a_Rv, & -!$OMP hP_Lv,hP_Rv,I_maxitt,khdt_epi_y,Tr_flux_3d, & -!$OMP Tr_adj_vert_L,Tr_adj_vert_R,Idt) & -!$OMP private(Tr_min_face,Tr_max_face,kLa,kLb,kRa,kRb, & -!$OMP Tr_La,Tr_Lb,Tr_Ra,Tr_Rb,Tr_av_L,wt_b,Tr_av_R, & -!$OMP h_L,h_R,Tr_flux,Tr_adj_vert,wt_a,vol) + !$OMP parallel do default(shared) private(Tr_min_face,Tr_max_face,kLa,kLb,kRa,kRb, & + !$OMP Tr_La,Tr_Lb,Tr_Ra,Tr_Rb,Tr_av_L,wt_b,Tr_av_R, & + !$OMP h_L,h_R,Tr_flux,Tr_adj_vert,wt_a,vol) do J=js-1,je ; do i=is,ie ; if (G%mask2dCv(i,J) > 0.0) then ! Determine the fluxes through the meridional faces. @@ -1370,7 +1404,11 @@ subroutine tracer_epipycnal_ML_diff(h, dt, Tr, ntr, khdt_epi_x, khdt_epi_y, G, & kRb = kRa ; if (max_kRho(i,j+1) < nz) kRb = max_kRho(i,j+1)+1 Tr_La = Tr_min_face ; Tr_Lb = Tr_La ; Tr_Ra = Tr_La ; Tr_Rb = Tr_La if (h(i,j,kLa) > h_exclude) Tr_La = Tr(m)%t(i,j,kLa) - if (h(i,j,kLb) > h_exclude) Tr_La = Tr(m)%t(i,j,kLb) + if ((CS%answer_date <= 20240330) .and. CS%limit_bug) then + if (h(i,j,kLb) > h_exclude) Tr_La = Tr(m)%t(i,j,kLb) + else + if (h(i,j,kLb) > h_exclude) Tr_Lb = Tr(m)%t(i,j,kLb) + endif if (h(i,j+1,kRa) > h_exclude) Tr_Ra = Tr(m)%t(i,j+1,kRa) if (h(i,j+1,kRb) > h_exclude) Tr_Rb = Tr(m)%t(i,j+1,kRb) Tr_min_face = min(Tr_min_face, Tr_La, Tr_Lb, Tr_Ra, Tr_Rb) @@ -1464,42 +1502,69 @@ subroutine tracer_epipycnal_ML_diff(h, dt, Tr, ntr, khdt_epi_x, khdt_epi_y, G, & Tr(m)%df2d_y(i,J) = Tr(m)%df2d_y(i,J) + Tr_flux * Idt enddo ! Loop over pairings at faces. endif ; enddo ; enddo ! i- & j- loops over meridional faces. -!$OMP parallel do default(none) shared(is,ie,js,je,G,nPv,k0b_Lv,k0b_Rv,deep_wt_Lv, & -!$OMP tr_flux_conv,Tr_flux_3d,k0a_Lv,Tr_adj_vert_L,& -!$OMP deep_wt_Rv,k0a_Rv,Tr_adj_vert_R) & -!$OMP private(kLa,kLb,kRa,kRb,wt_b,wt_a) - do i=is,ie ; do J=js-1,je ; if (G%mask2dCv(i,J) > 0.0) then + + !$OMP parallel do default(shared) private(kLa,kLb,kRa,kRb,wt_b,wt_a) + do J=js-1,je ; do i=is,ie ; if (G%mask2dCv(i,J) > 0.0) then ! The non-stride-1 loop order here is to facilitate openMP threading. However, it might be ! suboptimal when openMP threading is not used, at which point it might be better to fuse - ! these loope with those that precede it and thereby eliminate the need for three 3-d arrays. - do k=1,nPv(i,J) - kLb = k0b_Lv(J)%p(i,k); kRb = k0b_Rv(J)%p(i,k) - if (deep_wt_Lv(J)%p(i,k) >= 1.0) then - tr_flux_conv(i,j,kLb) = tr_flux_conv(i,j,kLb) - Tr_flux_3d(i,J,k) - else - kLa = k0a_Lv(J)%p(i,k) - wt_b = deep_wt_Lv(J)%p(i,k) ; wt_a = 1.0 - wt_b - tr_flux_conv(i,j,kLa) = tr_flux_conv(i,j,kLa) - (wt_a*Tr_flux_3d(i,J,k) + Tr_adj_vert_L(i,J,k)) - tr_flux_conv(i,j,kLb) = tr_flux_conv(i,j,kLb) - (wt_b*Tr_flux_3d(i,J,k) - Tr_adj_vert_L(i,J,k)) - endif - if (deep_wt_Rv(J)%p(i,k) >= 1.0) then - tr_flux_conv(i,j+1,kRb) = tr_flux_conv(i,j+1,kRb) + Tr_flux_3d(i,J,k) - else - kRa = k0a_Rv(J)%p(i,k) - wt_b = deep_wt_Rv(J)%p(i,k) ; wt_a = 1.0 - wt_b - tr_flux_conv(i,j+1,kRa) = tr_flux_conv(i,j+1,kRa) + & - (wt_a*Tr_flux_3d(i,J,k) - Tr_adj_vert_R(i,J,k)) - tr_flux_conv(i,j+1,kRb) = tr_flux_conv(i,j+1,kRb) + & - (wt_b*Tr_flux_3d(i,J,k) + Tr_adj_vert_R(i,J,k)) - endif - enddo + ! this loop with those that precede it and thereby eliminate the need for three 3-d arrays. + if (CS%answer_date <= 20240330) then + do k=1,nPv(i,J) + kLb = k0b_Lv(J)%p(i,k); kRb = k0b_Rv(J)%p(i,k) + if (deep_wt_Lv(J)%p(i,k) >= 1.0) then + tr_flux_conv(i,j,kLb) = tr_flux_conv(i,j,kLb) - Tr_flux_3d(i,J,k) + else + kLa = k0a_Lv(J)%p(i,k) + wt_b = deep_wt_Lv(J)%p(i,k) ; wt_a = 1.0 - wt_b + tr_flux_conv(i,j,kLa) = tr_flux_conv(i,j,kLa) - (wt_a*Tr_flux_3d(i,J,k) + Tr_adj_vert_L(i,J,k)) + tr_flux_conv(i,j,kLb) = tr_flux_conv(i,j,kLb) - (wt_b*Tr_flux_3d(i,J,k) - Tr_adj_vert_L(i,J,k)) + endif + if (deep_wt_Rv(J)%p(i,k) >= 1.0) then + tr_flux_conv(i,j+1,kRb) = tr_flux_conv(i,j+1,kRb) + Tr_flux_3d(i,J,k) + else + kRa = k0a_Rv(J)%p(i,k) + wt_b = deep_wt_Rv(J)%p(i,k) ; wt_a = 1.0 - wt_b + tr_flux_conv(i,j+1,kRa) = tr_flux_conv(i,j+1,kRa) + & + (wt_a*Tr_flux_3d(i,J,k) - Tr_adj_vert_R(i,J,k)) + tr_flux_conv(i,j+1,kRb) = tr_flux_conv(i,j+1,kRb) + & + (wt_b*Tr_flux_3d(i,J,k) + Tr_adj_vert_R(i,J,k)) + endif + enddo + else + do k=1,nPv(i,J) + kLb = k0b_Lv(J)%p(i,k); kRb = k0b_Rv(J)%p(i,k) + if (deep_wt_Lv(J)%p(i,k) >= 1.0) then + tr_flux_N(i,j,kLb) = tr_flux_N(i,j,kLb) + Tr_flux_3d(i,J,k) + else + kLa = k0a_Lv(J)%p(i,k) + wt_b = deep_wt_Lv(J)%p(i,k) ; wt_a = 1.0 - wt_b + tr_flux_N(i,j,kLa) = tr_flux_N(i,j,kLa) + (wt_a*Tr_flux_3d(i,J,k) + Tr_adj_vert_L(i,J,k)) + tr_flux_N(i,j,kLb) = tr_flux_N(i,j,kLb) + (wt_b*Tr_flux_3d(i,J,k) - Tr_adj_vert_L(i,J,k)) + endif + if (deep_wt_Rv(J)%p(i,k) >= 1.0) then + tr_flux_S(i,j+1,kRb) = tr_flux_S(i,j+1,kRb) + Tr_flux_3d(i,J,k) + else + kRa = k0a_Rv(J)%p(i,k) + wt_b = deep_wt_Rv(J)%p(i,k) ; wt_a = 1.0 - wt_b + tr_flux_S(i,j+1,kRa) = tr_flux_S(i,j+1,kRa) + (wt_a*Tr_flux_3d(i,J,k) - Tr_adj_vert_R(i,J,k)) + tr_flux_S(i,j+1,kRb) = tr_flux_S(i,j+1,kRb) + (wt_b*Tr_flux_3d(i,J,k) + Tr_adj_vert_R(i,J,k)) + endif + enddo + endif endif ; enddo ; enddo + + if (CS%answer_date >= 20240331) then + !$OMP parallel do default(shared) + do k=1,PEmax_kRho ; do j=js,je ; do i=is,ie + tr_flux_conv(i,j,k) = ((tr_flux_W(i,j,k) - tr_flux_E(i,j,k)) + & + (tr_flux_S(i,j,k) - tr_flux_N(i,j,k))) + enddo ; enddo ; enddo + endif + !$OMP parallel do default(shared) do k=1,PEmax_kRho ; do j=js,je ; do i=is,ie if ((G%mask2dT(i,j) > 0.0) .and. (h(i,j,k) > 0.0)) then - Tr(m)%t(i,j,k) = Tr(m)%t(i,j,k) + tr_flux_conv(i,j,k) / & - (h(i,j,k)*G%areaT(i,j)) - tr_flux_conv(i,j,k) = 0.0 + Tr(m)%t(i,j,k) = Tr(m)%t(i,j,k) + tr_flux_conv(i,j,k) / (h(i,j,k)*G%areaT(i,j)) endif enddo ; enddo ; enddo @@ -1546,6 +1611,7 @@ subroutine tracer_hor_diff_init(Time, G, GV, US, param_file, diag, EOS, diabatic ! This include declares and sets the variable "version". # include "version_variable.h" character(len=40) :: mdl = "MOM_tracer_hor_diff" ! This module's name. + integer :: default_answer_date if (associated(CS)) then call MOM_error(WARNING, "tracer_hor_diff_init called with associated control structure.") @@ -1604,6 +1670,21 @@ subroutine tracer_hor_diff_init(Time, G, GV, US, param_file, diag, EOS, diabatic "If true, then recalculate the neutral surfaces if the \n"//& "diffusive CFL is exceeded. If false, assume that the \n"//& "positions of the surfaces do not change \n", default=.false.) + call get_param(param_file, mdl, "DEFAULT_ANSWER_DATE", default_answer_date, & + "This sets the default value for the various _ANSWER_DATE parameters.", & + default=99991231, do_not_log=.true.) + call get_param(param_file, mdl, "HOR_DIFF_ANSWER_DATE", CS%answer_date, & + "The vintage of the order of arithmetic to use for the tracer diffusion. "//& + "Values of 20240330 or below recover the answers from the original form of the "//& + "along-isopycnal mixed layer to interior mixing code, while higher values use "//& + "mathematically equivalent expressions that recover rotational symmetry "//& + "when DIFFUSE_ML_TO_INTERIOR is true.", & + default=20240101, do_not_log=.not.CS%Diffuse_ML_interior) + !### Change the default later to default_answer_date. + call get_param(param_file, mdl, "HOR_DIFF_LIMIT_BUG", CS%limit_bug, & + "If true and the answer date is 20240330 or below, use a rotational symmetry "//& + "breaking bug when limiting the tracer properties in tracer_epipycnal_ML_diff.", & + default=.true., do_not_log=((.not.CS%Diffuse_ML_interior).or.(CS%answer_date>=20240331))) CS%ML_KhTR_scale = 1.0 if (CS%Diffuse_ML_interior) then call get_param(param_file, mdl, "ML_KHTR_SCALE", CS%ML_KhTR_scale, &