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

Changes redi tapering and slopes #65

Closed
wants to merge 4 commits into from
Closed

Changes redi tapering and slopes #65

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
0301b8e
Changes redi tapering and slopes
vanroekel Sep 6, 2023
51fe51f
fixes build errors
vanroekel Sep 6, 2023
a230f47
changes redi namelist options
vanroekel Sep 6, 2023
5eb06c8
Fixed OpenMP
mark-petersen Sep 14, 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
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
12 changes: 11 additions & 1 deletion 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 @@ -127,7 +127,7 @@ module ocn_diagnostics_variables
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 @@ -420,6 +420,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 +723,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 +734,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 @@ -961,6 +967,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 @@ -1161,6 +1169,8 @@ subroutine ocn_diagnostics_variables_destroy(err) !{{{
SSHGradient, &
circulation, &
slopeTriadUp, &
limiterUp, &
limiterDown, &
gradSSHZonal, &
cGMphaseSpeed, &
velocityZonal, &
Expand Down
82 changes: 40 additions & 42 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 @@ -229,30 +245,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 +310,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 +355,34 @@ 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
RediKappaScaling(1:minLevelCell(iCell),iCell) = 0.0_RKIND
RediKappaScaling(maxLevelCell(iCell):nVertLevels,iCell) = 0.0_RKIND
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 +900,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 +914,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
11 changes: 5 additions & 6 deletions components/mpas-ocean/src/shared/mpas_ocn_tracer_hmix.F
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -80,8 +80,8 @@ module ocn_tracer_hmix
!-----------------------------------------------------------------------

subroutine ocn_tracer_hmix_tend(meshPool, layerThickEdgeMean, layerThickness, zMid, tracers, &
RediKappa, slopeTriadUp, slopeTriadDown, dt, isActiveTracer, &
RediKappaScaling, rediLimiterCount, tend, err)!{{{
RediKappa, slopeTriadUp, slopeTriadDown, limiterUp, limiterDown, &
dt, isActiveTracer, rediLimiterCount, tend, err)!{{{


!-----------------------------------------------------------------
Expand All @@ -98,11 +98,9 @@ subroutine ocn_tracer_hmix_tend(meshPool, layerThickEdgeMean, layerThickness, zM
layerThickEdgeMean, &!< Input: mean thickness at edge
layerThickness, &!< Input: thickness at centers
zMid !< Input: Z coordinate at the center of a cell
real (kind=RKIND), dimension(:,:), intent(in), optional :: &
RediKappaScaling !< Input: vertical structure of GM/Redi limiter

real (kind=RKIND), dimension(:,:,:), intent(in) :: &
slopeTriadUp, slopeTriadDown, &
slopeTriadUp, slopeTriadDown, limiterUp, limiterDown, &
tracers !< Input: tracer quantities

real (kind=RKIND), intent(in) :: dt
Expand Down Expand Up @@ -154,7 +152,8 @@ subroutine ocn_tracer_hmix_tend(meshPool, layerThickEdgeMean, layerThickness, zM
call ocn_tracer_hmix_Redi_tend(meshPool, layerThickness, layerThickEdgeMean, zMid, tracers, &
RediKappa, dt, isActiveTracer, &
slopeTriadUp, slopeTriadDown, &
RediKappaScaling, rediLimiterCount, tend, err1)
limiterUp, limiterDown, rediLimiterCount, &
tend, err1)
endif
err = ior(err1, err2)
call mpas_timer_stop("tracer hmix")
Expand Down
Loading