diff --git a/src/tracer/MOM_tracer_hor_diff.F90 b/src/tracer/MOM_tracer_hor_diff.F90 index 6f4e5d0f90..920c03d484 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 20240309 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 20240309 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 20240309 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 20240309 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, &