diff --git a/cime_config/namelist_definition_blom.xml b/cime_config/namelist_definition_blom.xml
index 3dc4f527..dfd114e2 100644
--- a/cime_config/namelist_definition_blom.xml
+++ b/cime_config/namelist_definition_blom.xml
@@ -1853,6 +1853,17 @@
Type of lateral tracer eddy diffusion: Valid methods:
+
+ logical
+ diffusion
+ diffusion
+
+ .false.
+
+ If true, layers constructed for diffusive flux computations are gradually aligned with the surface within the mixed layer
+
+
+
diff --git a/cime_config/ocn_in.readme b/cime_config/ocn_in.readme
index 245f3751..15c65228 100644
--- a/cime_config/ocn_in.readme
+++ b/cime_config/ocn_in.readme
@@ -182,6 +182,9 @@
! types: 'none', 'vr12-ma', 'lf17'
! LTEDTP : Type of lateral tracer eddy diffusion: Valid methods:
! 'layer', 'neutral'.
+! NDIFF_SURFACE_ALIGN : If true, layers constructed for diffusive flux
+! computations are gradually aligned with the surface
+! within the mixed layer.
!===========================================================================
! NAMELIST FOR CHANNEL WIDTH MODIFICATIONS
diff --git a/phy/mod_ale_regrid_remap.F90 b/phy/mod_ale_regrid_remap.F90
index 23ff95a5..7d7d9f21 100644
--- a/phy/mod_ale_regrid_remap.F90
+++ b/phy/mod_ale_regrid_remap.F90
@@ -45,10 +45,11 @@ module mod_ale_regrid_remap
extract_polycoeff, regrid, &
prepare_remapping, remap, &
hor3map_noerr, hor3map_errstr
- use mod_diffusion, only: ltedtp_opt, ltedtp_neutral, difmxp, &
- utflld, vtflld, usflld, vsflld
+ use mod_diffusion, only: ltedtp_opt, ltedtp_neutral, ndiff_surface_align, &
+ difmxp, utflld, vtflld, usflld, vsflld
use mod_ndiff, only: ndiff_prep_jslice, ndiff_uflx_jslice, &
ndiff_vflx_jslice, ndiff_update_trc_jslice
+ use mod_cmnfld, only: dpml
use mod_dia, only: ddm, nphy, alarm_phy, &
depthslev_bnds, pbath, ubath, vbath, phylvl, &
acc_templvl, acc_salnlvl, &
@@ -1603,6 +1604,7 @@ subroutine ale_regrid_remap(m, n, mm, nn, k1m, k1n)
do nt = 1, ntr
call xctilr(trc(1-nbdy,1-nbdy,k1n,nt), 1, kk, 1, 1, halo_ps)
enddo
+ if (ndiff_surface_align) call xctilr(dpml, 1, 1, 1, 1, halo_ps)
end if
! Initial j-slice indices.
diff --git a/phy/mod_diffusion.F90 b/phy/mod_diffusion.F90
index e5408bd1..0148e7b4 100644
--- a/phy/mod_diffusion.F90
+++ b/phy/mod_diffusion.F90
@@ -77,8 +77,11 @@ module mod_diffusion
! diffusivities.
bdmldp, & ! If true, make the background mixing latitude dependent
! according to Gregg et al. (2003).
- smobld ! If true, apply lateral smoothing of CVMix estimated boundary
+ smobld, & ! If true, apply lateral smoothing of CVMix estimated boundary
! layer depth.
+ ndiff_surface_align ! If true, layers constructed for diffusive flux
+ ! computations are gradually aligned with the surface
+ ! within the mixed layer.
character(len = fnmlen) :: &
tbfile ! Name of file containing topographic beta parameter.
character(len = 80) :: &
@@ -180,9 +183,10 @@ module mod_diffusion
public :: egc, eggam, eglsmn, egmndf, egmxdf, egidfq, rhiscf, ri0, &
bdmc1, bdmc2, bdmldp, iwdflg, iwdfac, nubmin, tkepf, bdmtyp, &
eddf2d, edsprs, edanis, redi3d, rhsctp, tbfile, edfsmo, smobld, &
- lngmtp, eitmth_opt, eitmth_intdif, eitmth_gm, edritp_opt, &
- edritp_shear, edritp_large_scale, edwmth_opt, edwmth_smooth, &
- edwmth_step, ltedtp_opt, ltedtp_layer, ltedtp_neutral, &
+ ndiff_surface_align, lngmtp, eitmth_opt, eitmth_intdif, &
+ eitmth_gm, edritp_opt, edritp_shear, edritp_large_scale, &
+ edwmth_opt, edwmth_smooth, edwmth_step, ltedtp_opt, ltedtp_layer, &
+ ltedtp_neutral, &
difint, difiso, difdia, difmxp, difmxq, difwgt, &
umfltd, vmfltd, umflsm, vmflsm, utfltd, vtfltd, &
utflsm, vtflsm, utflld, vtflld, usfltd, vsfltd, &
@@ -209,7 +213,7 @@ subroutine readnml_diffusion
egc, eggam, eglsmn, egmndf, egmxdf, egidfq, rhiscf, ri0, &
bdmc1, bdmc2, bdmldp, iwdflg, iwdfac, nubmin, tkepf, bdmtyp, eddf2d, &
edsprs, edanis, redi3d, rhsctp, tbfile, edfsmo, smobld, lngmtp, &
- eitmth, edritp, edwmth, ltedtp
+ eitmth, edritp, edwmth, ltedtp, ndiff_surface_align
! Read variables in the namelist group 'diffusion'.
if (mnproc == 1) then
@@ -268,6 +272,7 @@ subroutine readnml_diffusion
call xcbcst(edritp)
call xcbcst(edwmth)
call xcbcst(ltedtp)
+ call xcbcst(ndiff_surface_align)
endif
if (mnproc == 1) then
write (lp,*) 'readnml_diffusion: diffusion variables:'
@@ -300,6 +305,7 @@ subroutine readnml_diffusion
write (lp,*) ' edritp = ', trim(edritp)
write (lp,*) ' edwmth = ', trim(edwmth)
write (lp,*) ' ltedtp = ', trim(ltedtp)
+ write (lp,*) ' ndiff_surface_align = ', ndiff_surface_align
endif
! Resolve options.
diff --git a/phy/mod_ndiff.F90 b/phy/mod_ndiff.F90
index 780a24b5..dff4a02a 100644
--- a/phy/mod_ndiff.F90
+++ b/phy/mod_ndiff.F90
@@ -29,8 +29,9 @@ module mod_ndiff
use mod_grid, only: scuy, scvx, scp2, scuxi, scvyi
use mod_eos, only: drhodt, drhods, rho
use mod_state, only: dp, temp, saln, utflx, vtflx, usflx, vsflx, pu, pv
- use mod_diffusion, only: difiso, utflld, vtflld, usflld, vsflld
- use mod_cmnfld, only: nslpx, nslpy
+ use mod_diffusion, only: ndiff_surface_align, difiso, utflld, vtflld, &
+ usflld, vsflld
+ use mod_cmnfld, only: nslpx, nslpy, dpml
use mod_tracers, only: trc
implicit none
@@ -192,18 +193,19 @@ pure subroutine ndiff_flx(p_srcdi_m, t_srcdi_m, tpc_src_m, &
real(r8), dimension(ntr_loc,2) :: t_ni_m, t_ni_p
real(r8), dimension(ntr_loc) :: t_nl_m, t_nl_p
real(r8), dimension(2) :: x_ni_m, x_ni_p, p_ni_m, p_ni_p
- real(r8) :: drho_curr, p_ni_m_prev, p_ni_p_prev, &
- drhodt_x0, drhodt_x1, drhods_x0, drhods_x1, &
- x_ni, p_ni, drho_prev, dp_dst_u, dp_dst_l, pu_m, pl_m, pu_p, pl_p, &
- pp1, pp2, dp_ni_m, dp_ni_p, dp_ni, q, dt, ds, &
- tflx, sflx, p_ni_up, p_ni_lo, dp_ni_i, mlfrac, p_nslp_dst
- integer :: nns, is_m, is_p, ks_m, ks_p, k, ks_m_prev, ks_p_prev, &
- kd_m, kd_p, isn_m, isn_p, ksn_m, ksn_p, &
- nip, nic, kuv, case_m, case_p, nt, kuvm, kd, ks
+ real(r8) :: pml, drho_curr, p_ni_m_prev, p_ni_p_prev, &
+ drhodt_x0, drhodt_x1, drhods_x0, drhods_x1, &
+ x_ni, p_ni, drho_prev, p1_m, p2_m, p1_p, p2_p, &
+ dp_dst_u, dp_dst_l, pu_m, pl_m, pu_p, pl_p, &
+ pp1, pp2, dp_ni_m, dp_ni_p, dp_ni, q, dt, ds, &
+ tflx, sflx, p_ni_up, p_ni_lo, dp_ni_i, mlfrac, p_nslp_dst
+ integer :: nns, issa_m, issa_p, kssa_m, kssa_p, is_m, is_p, ks_m, ks_p, k, &
+ ks_m_prev, ks_p_prev, kd_m, kd_p, isn_m, isn_p, ksn_m, ksn_p, &
+ nip, nic, kuv, case_m, case_p, nt, kuvm, kd, ks
logical, dimension(kk) :: stab_src_m, stab_src_p
logical :: drho_neg, drho_pos, drho_zero, &
- advance_src_m, advance_src_p, advance_dst_m, advance_dst_p, &
- found_ni
+ advance_src_m, advance_src_p, advance_dst_m, advance_dst_p, &
+ found_ni
real(r8), parameter :: mval = 1.e30_r8
@@ -221,22 +223,55 @@ pure subroutine ndiff_flx(p_srcdi_m, t_srcdi_m, tpc_src_m, &
nns = 0
- is_m = 1
- ks_m = 1
- ks_p = 1
- is_p = 1
- drho_curr = drho(t_srcdi_m(is_m,ks_m,it), &
- t_srcdi_m(is_m,ks_m,is), &
- t_srcdi_p(is_p,ks_p,it), &
- t_srcdi_p(is_p,ks_p,is), &
- .5_r8*( drhodt_srcdi_m(is_m,ks_m) &
- + drhodt_srcdi_p(is_p,ks_p)), &
- .5_r8*( drhods_srcdi_m(is_m,ks_m) &
- + drhods_srcdi_p(is_p,ks_p)))
- p_ni_m_prev = p_srcdi_m(1,1)
- p_ni_p_prev = p_srcdi_p(1,1)
-
- search_loop1: do
+ if (ndiff_surface_align) then
+
+ ! With ndiff_surface_align = .true., interfaces of layers constructed for
+ ! diffusive flux computations are gradually aligned with the surface
+ ! within the mixed layer. Therefore, start the search for neutral
+ ! interfaces for layer interfaces below the mixed layer depth.
+
+ pml = .5_r8*( p_srcdi_m(1,1) + dpml(i_m,j_m) &
+ + p_srcdi_p(1,1) + dpml(i_p,j_p))
+ kssa_m = 2
+ do while (kssa_m <= ksmx_m)
+ if (p_srcdi_m(1,kssa_m) > pml) exit
+ kssa_m = kssa_m + 1
+ enddo
+ kssa_p = 2
+ do while (kssa_p <= ksmx_p)
+ if (p_srcdi_p(1,kssa_p) > pml) exit
+ kssa_p = kssa_p + 1
+ enddo
+ is_m = 1
+ ks_m = kssa_m
+ is_p = 1
+ ks_p = kssa_p
+ p_ni_m_prev = pml
+ p_ni_p_prev = pml
+ else
+
+ ! Start the search for neutral interfaces from the uppermost source layer
+ ! interface.
+
+ is_m = 1
+ ks_m = 1
+ is_p = 1
+ ks_p = 1
+ p_ni_m_prev = p_srcdi_m(1,1)
+ p_ni_p_prev = p_srcdi_p(1,1)
+ endif
+
+ if (ks_m <= ksmx_m .and. ks_p <= ksmx_p) &
+ drho_curr = drho(t_srcdi_m(is_m,ks_m,it), &
+ t_srcdi_m(is_m,ks_m,is), &
+ t_srcdi_p(is_p,ks_p,it), &
+ t_srcdi_p(is_p,ks_p,is), &
+ .5_r8*( drhodt_srcdi_m(is_m,ks_m) &
+ + drhodt_srcdi_p(is_p,ks_p)), &
+ .5_r8*( drhods_srcdi_m(is_m,ks_m) &
+ + drhods_srcdi_p(is_p,ks_p)))
+
+ search_loop1: do while (ks_m <= ksmx_m .and. ks_p <= ksmx_p)
drho_neg = drho_curr <= - rho_eps
drho_pos = drho_curr >= rho_eps
@@ -252,7 +287,7 @@ pure subroutine ndiff_flx(p_srcdi_m, t_srcdi_m, tpc_src_m, &
drhods_x0 = .5_r8*( drhods_srcdi_m(1 ,ks_m) &
+ drhods_srcdi_p(is_p,ks_p))
drhods_x1 = .5_r8*( drhods_srcdi_m(2 ,ks_m) &
- + drhods_srcdi_p(is_p,ks_p))
+ + drhods_srcdi_p(is_p,ks_p))
x_ni = drhoroot(tpc_src_m(:,ks_m,it), tpc_src_m(:,ks_m,is), &
t_srcdi_p(is_p,ks_p,it), t_srcdi_p(is_p,ks_p,is), &
drhodt_x1, drhodt_x0, drhods_x1, drhods_x0)
@@ -358,6 +393,81 @@ pure subroutine ndiff_flx(p_srcdi_m, t_srcdi_m, tpc_src_m, &
enddo search_loop1
+ if (ndiff_surface_align) then
+
+ ! Above the uppermost identified neutral interface, gradually align
+ ! interfaces of layers for diffusive flux computations with the surface.
+
+ issa_m = 1
+ do while (kssa_m <= ksmx_m)
+ if (p_ni_srcdi_m(issa_m,kssa_m) /= mval) exit
+ if (issa_m == 1) then
+ issa_m = 2
+ else
+ kssa_m = kssa_m + 1
+ issa_m = 1
+ endif
+ enddo
+ issa_p = 1
+ do while (kssa_p <= ksmx_p)
+ if (p_ni_srcdi_p(issa_p,kssa_p) /= mval) exit
+ if (issa_p == 1) then
+ issa_p = 2
+ else
+ kssa_p = kssa_p + 1
+ issa_p = 1
+ endif
+ enddo
+ if (kssa_m > ksmx_m .or. kssa_p > ksmx_p) then
+ p_ni_srcdi_m(1,1) = p_srcdi_m(1,1)
+ do ks_m = 1, ksmx_m - 1
+ if (p_srcdi_m(1,ks_m) > p_srcdi_p(2,ksmx_p)) exit
+ p_ni = min(p_srcdi_m(2,ks_m), p_srcdi_p(2,ksmx_p))
+ p_ni_srcdi_m(1,ks_m+1) = p_ni
+ p_ni_srcdi_m(2,ks_m ) = p_ni
+ stab_src_m(ks_m) = .true.
+ enddo
+ p_ni_srcdi_p(1,1) = p_srcdi_p(1,1)
+ do ks_p = 1, ksmx_p - 1
+ if (p_srcdi_p(1,ks_p) > p_srcdi_m(2,ksmx_m)) exit
+ p_ni = min(p_srcdi_p(2,ks_p), p_srcdi_m(2,ksmx_m))
+ p_ni_srcdi_p(1,ks_p+1) = p_ni
+ p_ni_srcdi_p(2,ks_p ) = p_ni
+ stab_src_p(ks_p) = .true.
+ enddo
+ else
+ if (p_srcdi_m(issa_m,kssa_m) < p_ni_srcdi_p(issa_p,kssa_p)) then
+ p1_m = p_srcdi_m(1,1)
+ p2_m = p_srcdi_m(issa_m,kssa_m)
+ p1_p = p_srcdi_p(1,1)
+ p2_p = p_ni_srcdi_m(issa_m,kssa_m)
+ else
+ p1_m = p_srcdi_m(1,1)
+ p2_m = p_ni_srcdi_p(issa_p,kssa_p)
+ p1_p = p_srcdi_p(1,1)
+ p2_p = p_srcdi_p(issa_p,kssa_p)
+ endif
+ p_ni_srcdi_m(1,1) = p1_p
+ do ks_m = 1, kssa_m - 1
+ p_ni = ( (p_srcdi_m(2,ks_m) - p1_m)*p2_p &
+ + (p2_m - p_srcdi_m(2,ks_m))*p1_p) &
+ /(p2_m - p1_m)
+ p_ni_srcdi_m(1,ks_m+1) = p_ni
+ p_ni_srcdi_m(2,ks_m ) = p_ni
+ stab_src_m(ks_m) = .true.
+ enddo
+ p_ni_srcdi_p(1,1) = p1_m
+ do ks_p = 1, kssa_p -1
+ p_ni = ( (p_srcdi_p(2,ks_p) - p1_p)*p2_m &
+ + (p2_p - p_srcdi_p(2,ks_p))*p1_m) &
+ /(p2_p - p1_p)
+ p_ni_srcdi_p(1,ks_p+1) = p_ni
+ p_ni_srcdi_p(2,ks_p ) = p_ni
+ stab_src_p(ks_p) = .true.
+ enddo
+ endif
+ endif
+
! ------------------------------------------------------------------------
! Do another search from the surface, this time including destination
! interfaces, to identify neutral layers and compute fluxes that are added
@@ -434,8 +544,7 @@ pure subroutine ndiff_flx(p_srcdi_m, t_srcdi_m, tpc_src_m, &
ks_m = ks_m + 1
if (ks_m > ksmx_m) exit search_loop2
is_m = 1
- if (stab_src_m(ks_m) .and. p_ni_srcdi_m(is_m,ks_m) /= mval) &
- exit
+ if (stab_src_m(ks_m) .and. p_ni_srcdi_m(is_m,ks_m) /= mval) exit
endif
enddo
isn_m = is_m
@@ -460,9 +569,7 @@ pure subroutine ndiff_flx(p_srcdi_m, t_srcdi_m, tpc_src_m, &
ks_p = ks_p + 1
if (ks_p > ksmx_p) exit search_loop2
is_p = 1
- if (stab_src_p(ks_p) .and. p_ni_srcdi_p(is_p,ks_p) /= mval) then
- exit
- end if
+ if (stab_src_p(ks_p) .and. p_ni_srcdi_p(is_p,ks_p) /= mval) exit
endif
enddo
isn_p = is_p
@@ -503,14 +610,12 @@ pure subroutine ndiff_flx(p_srcdi_m, t_srcdi_m, tpc_src_m, &
if (kd_p > kdmx_p) exit search_loop2
endif
- do while (p_dstsnp_m(kd_m+1) &
- <= max(p_srcdi_m(1,ks_m), p_ni_m(nip)))
+ do while (p_dstsnp_m(kd_m+1) <= max(p_srcdi_m(1,ks_m), p_ni_m(nip)))
kd_m = kd_m + 1
if (kd_m > kdmx_m) exit search_loop2
enddo
- do while (p_dstsnp_p(kd_p+1) &
- <= max(p_srcdi_p(1,ks_p), p_ni_p(nip)))
+ do while (p_dstsnp_p(kd_p+1) <= max(p_srcdi_p(1,ks_p), p_ni_p(nip)))
kd_p = kd_p + 1
if (kd_p > kdmx_p) exit search_loop2
enddo
diff --git a/tests/fuk95/limits b/tests/fuk95/limits
index d11274b4..7ccb3160 100644
--- a/tests/fuk95/limits
+++ b/tests/fuk95/limits
@@ -300,6 +300,10 @@
! types: 'none', 'vr12-ma', 'lf17'
! LTEDTP : Type of lateral tracer eddy diffusion: Valid methods:
! 'layer', 'neutral'.
+! NDIFF_SURFACE_ALIGN : If true, layers constructed for diffusive flux
+! computations are gradually aligned with the surface
+! within the mixed layer.
+
&DIFFUSION
EITMTH = 'gm'
EDRITP = 'large scale'
@@ -330,6 +334,7 @@
SMOBLD = .true.
LNGMTP = 'none'
LTEDTP = 'layer'
+ NDIFF_SURFACE_ALIGN = .false.
/
! NAMELIST FOR CHANNEL WIDTH MODIFICATIONS