Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Improve slope calculation and fix slope limiting for redi mixing #5947

Merged
merged 6 commits into from
Oct 26, 2023
Merged

Improve slope calculation and fix slope limiting for redi mixing #5947

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
6 commits
Select commit Hold shift + click to select a range
b031624
Changes redi tapering and slopes
vanroekel Sep 6, 2023
135dea2
fixes build errors
vanroekel Sep 6, 2023
fbb9ce2
changes redi namelist options
vanroekel Sep 6, 2023
a0fdeca
Fixed OpenMP
mark-petersen Sep 14, 2023
fac4c9c
Update bld files to match Registry changes
jonbob Sep 22, 2023
534209c
Remove RediKappaScaling creation and initialization
mark-petersen Sep 25, 2023
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
4 changes: 1 addition & 3 deletions components/mpas-ocean/bld/build-namelist
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -557,9 +557,7 @@ add_default($nl, 'config_Redi_constant_kappa');
add_default($nl, 'config_Redi_maximum_slope');
add_default($nl, 'config_Redi_use_slope_taper');
add_default($nl, 'config_Redi_use_surface_taper');
add_default($nl, 'config_Redi_N2_based_taper_enable');
add_default($nl, 'config_Redi_N2_based_taper_min');
add_default($nl, 'config_Redi_N2_based_taper_limit_term1');
add_default($nl, 'config_Redi_limit_term1');
add_default($nl, 'config_Redi_min_layers_diag_terms');
add_default($nl, 'config_Redi_horizontal_taper');
add_default($nl, 'config_Redi_horizontal_ramp_min');
Expand Down
4 changes: 1 addition & 3 deletions components/mpas-ocean/bld/build-namelist-section
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -96,9 +96,7 @@ add_default($nl, 'config_Redi_constant_kappa');
add_default($nl, 'config_Redi_maximum_slope');
add_default($nl, 'config_Redi_use_slope_taper');
add_default($nl, 'config_Redi_use_surface_taper');
add_default($nl, 'config_Redi_N2_based_taper_enable');
add_default($nl, 'config_Redi_N2_based_taper_min');
add_default($nl, 'config_Redi_N2_based_taper_limit_term1');
add_default($nl, 'config_Redi_limit_term1');
add_default($nl, 'config_Redi_min_layers_diag_terms');
add_default($nl, 'config_Redi_horizontal_taper');
add_default($nl, 'config_Redi_horizontal_ramp_min');
Expand Down
4 changes: 1 addition & 3 deletions components/mpas-ocean/bld/namelist_files/namelist_defaults_mpaso.xml
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -141,9 +141,7 @@
<config_Redi_maximum_slope>0.01</config_Redi_maximum_slope>
<config_Redi_use_slope_taper>.true.</config_Redi_use_slope_taper>
<config_Redi_use_surface_taper>.true.</config_Redi_use_surface_taper>
<config_Redi_N2_based_taper_enable>.true.</config_Redi_N2_based_taper_enable>
<config_Redi_N2_based_taper_min>0.1</config_Redi_N2_based_taper_min>
<config_Redi_N2_based_taper_limit_term1>.true.</config_Redi_N2_based_taper_limit_term1>
<config_Redi_limit_term1>.true.</config_Redi_limit_term1>
<config_Redi_min_layers_diag_terms>6</config_Redi_min_layers_diag_terms>
<config_Redi_min_layers_diag_terms ocn_grid="ARRM10to60E2r1">15</config_Redi_min_layers_diag_terms>
<config_Redi_horizontal_taper>'ramp'</config_Redi_horizontal_taper>
Expand Down
18 changes: 1 addition & 17 deletions components/mpas-ocean/bld/namelist_files/namelist_definition_mpaso.xml
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -423,23 +423,7 @@ Valid values: .true. or .false.
Default: Defined in namelist_defaults.xml
</entry>

<entry id="config_Redi_N2_based_taper_enable" type="logical"
category="Redi_isopycnal_mixing" group="Redi_isopycnal_mixing">
If true apply the N2 limiting of Danabasoglu and Marshall 2007

Valid values: .true. or .false.
Default: Defined in namelist_defaults.xml
</entry>

<entry id="config_Redi_N2_based_taper_min" type="real"
category="Redi_isopycnal_mixing" group="Redi_isopycnal_mixing">
minimum taper value for the N2 limiting of Danabasoglu and Marshall 2007

Valid values: greater than zero and less than 1
Default: Defined in namelist_defaults.xml
</entry>

<entry id="config_Redi_N2_based_taper_limit_term1" type="logical"
<entry id="config_Redi_limit_term1" type="logical"
category="Redi_isopycnal_mixing" group="Redi_isopycnal_mixing">
If true, the N2 limiting is applied to the horizontal diffusion term

Expand Down
22 changes: 9 additions & 13 deletions components/mpas-ocean/src/Registry.xml
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -314,15 +314,7 @@
description="If true, Redi slope is tapered near sfc based on Large et al 1997"
possible_values=".true. or .false."
/>
<nml_option name="config_Redi_N2_based_taper_enable" type="logical" default_value=".true."
description="If true apply the N2 limiting of Danabasoglu and Marshall 2007"
possible_values=".true. or .false."
/>
<nml_option name="config_Redi_N2_based_taper_min" type="real" default_value="0.1"
description="minimum taper value for the N2 limiting of Danabasoglu and Marshall 2007"
possible_values="greater than zero and less than 1"
/>
<nml_option name="config_Redi_N2_based_taper_limit_term1" type="logical" default_value=".true."
<nml_option name="config_Redi_limit_term1" type="logical" default_value=".true."
description="If true, the N2 limiting is applied to the horizontal diffusion term"
possible_values=".true. or .false."
/>
Expand Down Expand Up @@ -3008,10 +3000,6 @@
description="spatially and depth varying GM kappa. The scaling is based on the Brunt Vaisala Frequency relative to a maximum value below the mixed layer, follows from Danabasoglu and Marshall 2007. If config_GM_closure is not set to N2_dependent the scaling value is set to 1 everywhere."
packages="gm"
/>
<var name="RediKappaScaling" type="real" dimensions="nVertLevelsP1 nCells Time" units="1"
description="Scaling coefficient for GM kappa. Varies from 0 to 1. The scaling is based on the Brunt Vaisala Frequency relative to a maximum value below the mixed layer, follows from Danabasoglu and Marshall 2007. If config_Redi_N2_based_taper_enable = .false. the scaling is set to 1 everywhere."
packages="gm"
/>
<var name="RediKappaSfcTaper" type="real" dimensions="nVertLevelsP1 nCells Time" units="1"
description="Scaling coefficient for Redi kappa. Varies from 0 to 1."
packages="gm"
Expand Down Expand Up @@ -3151,6 +3139,14 @@
description="Magnitude of slope of isopycnal surface, using triad through this cell and edge, angled up. Uses expansion of equation of state."
packages="forwardMode;analysisMode"
/>
<var name="limiterUp" type="real" dimensions="nVertLevels TWO nEdges Time" units="non-dimensional" default_value="0.0"
description="Magnitude of slope of isopycnal surface, using triad through this cell and edge, angled up. Uses expansion of equation of state."
packages="forwardMode;analysisMode"
/>
<var name="limiterDown" type="real" dimensions="nVertLevels TWO nEdges Time" units="non-dimensional" default_value="0.0"
description="Magnitude of slope of isopycnal surface, using triad through this cell and edge, angled up. Uses expansion of equation of state."
packages="forwardMode;analysisMode"
/>
<var name="k33" type="real" dimensions="nVertLevelsP1 nCells Time" units="m^2 s^-1"
description="The (3,3) entry of the Redi diffusion tensor. Added to the model vertical diffusion. Defined at the top of cell k"
packages="forwardMode;analysisMode"
Expand Down
18 changes: 11 additions & 7 deletions components/mpas-ocean/src/shared/mpas_ocn_diagnostics_variables.F
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -121,13 +121,12 @@ module ocn_diagnostics_variables
real (kind=RKIND), dimension(:,:), pointer :: velocityZonal

real(kind=RKIND), dimension(:,:), pointer :: gmKappaScaling
real(kind=RKIND), dimension(:,:), pointer :: RediKappaScaling
real(kind=RKIND), dimension(:,:), pointer :: RediKappaSfcTaper
real(kind=RKIND), dimension(:,:), pointer :: k33
real(kind=RKIND), dimension(:,:), pointer :: gradDensityEddy
real(kind=RKIND), dimension(:), pointer :: gmBolusKappa
real(kind=RKIND), dimension(:), pointer :: cGMphaseSpeed
real(kind=RKIND), dimension(:,:,:), pointer :: slopeTriadUp, slopeTriadDown
real(kind=RKIND), dimension(:,:,:), pointer :: limiterUp, limiterDown, slopeTriadUp, slopeTriadDown

integer, dimension(:), pointer :: indMLD
real(kind=RKIND), dimension(:), pointer :: dThreshMLD
Expand Down Expand Up @@ -265,8 +264,6 @@ subroutine ocn_diagnostics_variables_init(domain, jenkinsOn, hollandJenkinsOn, e
if (config_use_GM.or.config_use_Redi) then
call mpas_pool_get_array(diagnosticsPool, 'RediKappa', &
RediKappa)
call mpas_pool_get_array(diagnosticsPool, 'RediKappaScaling', &
RediKappaScaling)
call mpas_pool_get_array(diagnosticsPool, 'RediKappaSfcTaper', &
RediKappaSfcTaper)
call mpas_pool_get_array(diagnosticsPool, 'rediLimiterCount', &
Expand Down Expand Up @@ -420,6 +417,10 @@ subroutine ocn_diagnostics_variables_init(domain, jenkinsOn, hollandJenkinsOn, e
slopeTriadUp)
call mpas_pool_get_array(diagnosticsPool, 'slopeTriadDown', &
slopeTriadDown)
call mpas_pool_get_array(diagnosticsPool, 'limiterUp', &
limiterUp)
call mpas_pool_get_array(diagnosticsPool, 'limiterDown', &
limiterDown)

! GM-related fields
call mpas_pool_get_array(diagnosticsPool, 'cGMphaseSpeed', &
Expand Down Expand Up @@ -719,6 +720,7 @@ subroutine ocn_diagnostics_variables_init(domain, jenkinsOn, hollandJenkinsOn, e
!$acc SSHGradient, &
!$acc circulation, &
!$acc slopeTriadUp, &
!$acc limiterUp, &
!$acc gradSSHZonal, &
!$acc cGMphaseSpeed, &
!$acc velocityZonal, &
Expand All @@ -729,6 +731,7 @@ subroutine ocn_diagnostics_variables_init(domain, jenkinsOn, hollandJenkinsOn, e
!$acc layerThickEdgeFlux, &
!$acc layerThickEdgeMean, &
!$acc slopeTriadDown, &
!$acc limiterDown, &
!$acc normRelVortEdge, &
!$acc surfaceVelocity, &
!$acc vertVelocityTop, &
Expand Down Expand Up @@ -792,7 +795,6 @@ subroutine ocn_diagnostics_variables_init(domain, jenkinsOn, hollandJenkinsOn, e
!$acc RediKappa, &
!$acc gmBolusKappa, &
!$acc gmKappaScaling, &
!$acc RediKappaScaling, &
!$acc rediLimiterCount, &
!$acc RediKappaSfcTaper, &
!$acc RediHorizontalTaper, &
Expand Down Expand Up @@ -961,6 +963,8 @@ subroutine ocn_diagnostics_variables_destroy(err) !{{{
!$acc SSHGradient, &
!$acc circulation, &
!$acc slopeTriadUp, &
!$acc limiterUp, &
!$acc limiterDown, &
!$acc gradSSHZonal, &
!$acc cGMphaseSpeed, &
!$acc velocityZonal, &
Expand Down Expand Up @@ -1034,7 +1038,6 @@ subroutine ocn_diagnostics_variables_destroy(err) !{{{
!$acc RediKappa, &
!$acc gmBolusKappa, &
!$acc gmKappaScaling, &
!$acc RediKappaScaling, &
!$acc rediLimiterCount, &
!$acc RediKappaSfcTaper, &
!$acc RediHorizontalTaper, &
Expand Down Expand Up @@ -1161,6 +1164,8 @@ subroutine ocn_diagnostics_variables_destroy(err) !{{{
SSHGradient, &
circulation, &
slopeTriadUp, &
limiterUp, &
limiterDown, &
gradSSHZonal, &
cGMphaseSpeed, &
velocityZonal, &
Expand Down Expand Up @@ -1230,7 +1235,6 @@ subroutine ocn_diagnostics_variables_destroy(err) !{{{
nullify(RediKappa, &
gmBolusKappa, &
gmKappaScaling, &
RediKappaScaling, &
rediLimiterCount, &
RediKappaSfcTaper, &
RediHorizontalTaper, &
Expand Down
81 changes: 38 additions & 43 deletions components/mpas-ocean/src/shared/mpas_ocn_gm.F
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -62,6 +62,22 @@ module ocn_gm

contains

function safe_slope(num, den) result(slope)

real(kind=RKIND), intent(in) :: num, den
real(kind=RKIND) :: slope

! compute isopycnal slopes safely wrt. finite precision
! small slope assumption, so don't want |s| > 1

if (abs(den) > abs(num)) then
slope = num / den
else
slope = sign(1.0_RKIND, den * num)
end if

end function

!***********************************************************************
!
! routine ocn_GM_compute_Bolus_velocity
Expand Down Expand Up @@ -214,7 +230,6 @@ subroutine ocn_GM_compute_Bolus_velocity(statePool, &

!$omp do schedule(runtime)
do iCell = 1, nCells
RediKappaScaling(:, iCell) = 1.0_RKIND
RediKappaSfcTaper(:, iCell) = 1.0_RKIND
end do
!$omp end do
Expand All @@ -229,30 +244,6 @@ subroutine ocn_GM_compute_Bolus_velocity(statePool, &
allocate(dSdzTop(nVertLevels + 1))
allocate(k33Norm(nVertLevels + 1))

if ( config_Redi_N2_based_taper_enable ) then

!$omp parallel
!$omp do schedule(runtime) private(maxLocation, k, BruntVaisalaFreqTopEdge, maxN)
do iCell = 1, nCells
k = min(maxLevelCell(iCell) - 1, max(minLevelCell(iCell), indMLD(iCell)))
maxN = max(BruntVaisalaFreqTop(k, iCell), 0.0_RKIND)
do while (BruntVaisalaFreqTop(k + 1, iCell) > maxN .and. k < maxLevelCell(iCell) - 1)
k = k + 1
maxN = max(maxN,max(BruntVaisalaFreqTop(k, iCell), 0.0_RKIND))
enddo

maxLocation = k
do k = maxLocation, maxLevelCell(iCell)
BruntVaisalaFreqTopEdge = max(BruntVaisalaFreqTop(k, iCell), 0.0_RKIND)
RediKappaScaling(k, iCell) = min(max(config_Redi_N2_based_taper_min, &
BruntVaisalaFreqTopEdge/(maxN + 1.0E-10_RKIND)), &
1.0_RKIND)
end do
end do
!$omp end do
!$omp end parallel
end if

if (config_Redi_use_surface_taper) then
!$omp parallel
!$omp do schedule (runtime) private(zMLD, k)
Expand Down Expand Up @@ -318,17 +309,17 @@ subroutine ocn_GM_compute_Bolus_velocity(statePool, &
*dcEdgeInv
drhoDx = drhoDT*dTdx + drhoDS*dSdx


sfcTaperUp = drhoDT*dTdzTop(k) + drhoDS*dSdzTop(k)
sfcTaperDown = drhoDT*dTdzTop(k+1) + drhoDS*dSdzTop(k+1)

! Always compute *Up on the top cell and *Down on the bottom
! cell, even though they are never used. This avoids an if
! statement or separate computation for top and bottom.
slopeTriadUp(k, iCellSelf, iEdge) = &
-drhoDx/ &
(drhoDT*dTdzTop(k) &
+ drhoDS*dSdzTop(k) + 1E-15_RKIND)
safe_slope( -drhoDx, sfcTaperUp )
slopeTriadDown(k, iCellSelf, iEdge) = &
-drhoDx/ &
(drhoDT*dTdzTop(k + 1) &
+ drhoDS*dSdzTop(k + 1) + 1E-15_RKIND)
safe_slope( -drhoDx, sfcTaperDown )

! set taper of slope ('F' function from Danabasoglu and McWilliams 95)
if (abs(slopeTriadDown(k, iCellSelf, iEdge)) > 0.6_RKIND*config_redi_maximum_slope) then
Expand Down Expand Up @@ -363,30 +354,32 @@ subroutine ocn_GM_compute_Bolus_velocity(statePool, &
sfcTaperUp = 1.0_RKIND + sfcTaperFactor*(sfcTaper - 1.0_RKIND)
sfcTaperDown = 1.0_RKIND + sfcTaperFactor*(sfcTaper - 1.0_RKIND)

slopeTriadUp(k, iCellSelf, iEdge) = &
slopeTaperUp*sfcTaperUp*slopeTriadUp(k, iCellSelf, iEdge)
slopeTriadDown(k, iCellSelf, iEdge) = &
slopeTaperDown*sfcTaperDown*slopeTriadDown(k, iCellSelf, iEdge)

! Griffies 1998 eqn 34
! ! Griffies 1998 eqn 34
if (k > minLevelCell(iCell)) then
k33(k, iCell) = k33(k, iCell) + &
areaEdge*dzTop(k)*slopeTriadUp(k, iCellSelf, iEdge)**2
slopeTaperUp*sfcTaperUp*areaEdge*dzTop(k)*slopeTriadUp(k, iCellSelf, iEdge)**2
k33Norm(k) = k33Norm(k) + areaEdge*dzTop(k)
end if

!Taper Redi by tapering the slopes
! !Taper Redi by tapering the slopes
k33(k + 1, iCell) = k33(k + 1, iCell) + &
areaEdge*dzTop(k + 1)*slopeTriadDown(k, iCellSelf, iEdge)**2
slopeTaperDown*sfcTaperDown*areaEdge*dzTop(k + 1)*slopeTriadDown(k, iCellSelf, iEdge)**2
k33Norm(k + 1) = k33Norm(k + 1) + areaEdge*dzTop(k + 1)

limiterUp(k,iCellSelf,iEdge) = slopeTaperUp*sfcTaperUp
limiterDown(k,iCellSelf,iEdge) = slopeTaperDown*sfcTaperDown

slopeTriadUp(k, iCellSelf, iEdge) = &
slopeTaperUp*sfcTaperUp*slopeTriadUp(k, iCellSelf, iEdge)
slopeTriadDown(k, iCellSelf, iEdge) = &
slopeTaperDown*sfcTaperDown*slopeTriadDown(k, iCellSelf, iEdge)

end do ! maxLevelEdgeTop(iEdge)
end do ! nEdgesOnCell(iCell)

! Normalize k33
! ! Normalize k33
do k = minLevelCell(iCell)+1, maxLevelCell(iCell)
k33(k, iCell) = k33(k, iCell)/k33Norm(k)*RediKappaScaling(k, iCell)
k33(k, iCell) = k33(k, iCell)/k33Norm(k)
end do
k33(1:minLevelCell(iCell), iCell) = 0.0_RKIND
k33(maxLevelCell(iCell) + 1, iCell) = 0.0_RKIND
Expand Down Expand Up @@ -904,7 +897,7 @@ subroutine ocn_GM_init(domain, err)!{{{
real(kind=RKIND), dimension(:), pointer :: dcEdge

integer :: iEdge
integer, pointer :: nEdges
integer, pointer :: nCells, nEdges
integer, pointer :: nVertLevels
real(kind=RKIND), pointer :: sphere_radius
real(kind=RKIND), dimension(:), pointer :: latEdge, fEdge
Expand All @@ -918,10 +911,12 @@ subroutine ocn_GM_init(domain, err)!{{{
do while (associated(block))
call mpas_pool_get_subpool(block%structs, 'mesh', meshPool)
call mpas_pool_get_dimension(meshPool, 'nVertLevels', nVertLevels)
call mpas_pool_get_dimension(meshPool, 'nCells', nCells)
call mpas_pool_get_dimension(meshPool, 'nEdges', nEdges)
call mpas_pool_get_array(meshPool, 'cellsOnEdge', cellsOnEdge)
call mpas_pool_get_array(meshPool, 'dcEdge', dcEdge)
call mpas_pool_get_array(meshPool, 'fEdge', fEdge)

if (config_Redi_use_slope_taper) then
slopeTaperFactor = 1.0_RKIND
else
Expand Down
5 changes: 3 additions & 2 deletions components/mpas-ocean/src/shared/mpas_ocn_tendency.F
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1141,8 +1141,9 @@ subroutine ocn_tend_tracer(tendPool, statePool, forcingPool, &
call ocn_tracer_hmix_tend(meshPool, layerThickEdgeMean, &
layerThickness, zMid, tracerGroup, &
RediKappa, slopeTriadUp, &
slopeTriadDown, dt, isActiveTracer, &
RediKappaScaling, rediLimiterCount, &
slopeTriadDown, limiterUp, &
limiterDown, dt, isActiveTracer, &
rediLimiterCount, &
tracerGroupTend, err)

if (computeBudgets) then
Expand Down
Loading