Skip to content

Commit

Permalink
+Rotationally symmetric epipycnal diffusion option
Browse files Browse the repository at this point in the history 
  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.
  • Loading branch information
@Hallberg-NOAA @marshallward
Hallberg-NOAA authored and marshallward committed Apr 5, 2024
1 parent f9372f3 commit a3fd1f3
Showing 1 changed file with 149 additions and 68 deletions.
217 changes: 149 additions & 68 deletions src/tracer/MOM_tracer_hor_diff.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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.
Expand Down Expand Up @@ -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)
Expand Down Expand Up @@ -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]
Expand Down Expand Up @@ -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)
Expand Down Expand Up @@ -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

Expand All @@ -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.

Expand All @@ -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)
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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.

Expand All @@ -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)
Expand Down Expand Up @@ -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

Expand Down Expand Up @@ -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.")
Expand Down Expand Up @@ -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, &
Expand Down

0 comments on commit a3fd1f3

Please sign in to comment.