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

Fixes Redi Surface Taper double counting #5171

Merged
merged 3 commits into from
Sep 29, 2022
Merged

Fixes Redi Surface Taper double counting #5171

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
b547d67
Fixes Redi Surface Taper double counting
vanroekel Sep 2, 2022
98296b4
Removes extra characters
vanroekel Sep 2, 2022
42eca04
Add mixed layer depth check to surface taper
mark-petersen Sep 20, 2022
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
5 changes: 5 additions & 0 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 @@ -351,6 +351,11 @@ subroutine ocn_GM_compute_Bolus_velocity(statePool, &
slopeTaperDown = 1.0_RKIND + slopeTaperFactor*(slopeTaperDown - 1.0_RKIND)

sfcTaper = min(RediKappaSfcTaper(k, cell1), RediKappaSfcTaper(k, cell2))
if (k < min( indMLD(cell1), indMLD(cell2))) then
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
if (k < min( indMLD(cell1), indMLD(cell2))) then
if (k < min(indMLD(cell1), indMLD(cell2))) then

sfcTaper = 0.0_RKIND
else
sfcTaper = 1.0_RKIND
end if
Comment on lines 353 to +358
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I don't understand why the first line here isn't replaced by the last 5. It seems like the work done in the first line now does nothing.

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Also, could something be added to explain this change?

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

In particular, I'm wanting to make sure these still relate to the topic of this PR, the double counting of the Redi surface taper.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Your last point is a very good one. This is definitely beyond the title of the PR. This is now more along the lines of 'generic redi fixes'. I'll add a bit of explanation here about this change -- it is physically motivated and I believe fixes the stability issues @jonbob was seeing. Within the upper ocean the mesoscale eddies act horizontally and transition to along isopycnal in the interior (e.g. Ferrari et al 2005), this is a very rudimentary way to achieve that. The reason I think this will fix the issue with stability is that when I removed the surface taper double counting the redi cross terms became stronger, leading to numerical instability. This change is effectively a much stronger taper on those terms. I included this in my long 800 year spin up and didn't see stability issues, but Jon sees them quickly. I can add more explanation to this PR and/or in the code.

Finally, you are absolutely right that the first line of five is not necessary.

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Thanks! If you update the PR description, I am happy with this being included.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

more than happy to add this to the description, but I may hold off until @jonbob can test the impact to see if this stabilizes the run.


sfcTaperUp = 1.0_RKIND + sfcTaperFactor*(sfcTaper - 1.0_RKIND)
sfcTaperDown = 1.0_RKIND + sfcTaperFactor*(sfcTaper - 1.0_RKIND)
Expand Down
4 changes: 2 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 @@ -1136,8 +1136,8 @@ subroutine ocn_tend_tracer(tendPool, statePool, forcingPool, &
layerThickness, zMid, tracerGroup, &
RediKappa, slopeTriadUp, &
slopeTriadDown, dt, isActiveTracer, &
RediKappaSfcTaper, RediKappaScaling,&
rediLimiterCount, tracerGroupTend, err)
RediKappaScaling, rediLimiterCount, &
tracerGroupTend, err)

if (computeBudgets) then
! Finish budget calculation by subtracting before, after
Expand Down
6 changes: 3 additions & 3 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 @@ -81,7 +81,7 @@ module ocn_tracer_hmix

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


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

real (kind=RKIND), dimension(:,:,:), intent(in) :: &
slopeTriadUp, slopeTriadDown, &
Expand Down Expand Up @@ -153,7 +153,7 @@ subroutine ocn_tracer_hmix_tend(meshPool, layerThickEdgeMean, layerThickness, zM
if (config_use_Redi) then
call ocn_tracer_hmix_Redi_tend(meshPool, layerThickness, layerThickEdgeMean, zMid, tracers, &
RediKappa, dt, isActiveTracer, &
slopeTriadUp, slopeTriadDown, RediKappaSfcTaper, &
slopeTriadUp, slopeTriadDown, &
RediKappaScaling, rediLimiterCount, tend, err1)
endif
err = ior(err1, err2)
Expand Down
29 changes: 10 additions & 19 deletions components/mpas-ocean/src/shared/mpas_ocn_tracer_hmix_redi.F
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -75,7 +75,7 @@ module ocn_tracer_hmix_Redi

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

!-----------------------------------------------------------------
!
Expand All @@ -91,7 +91,6 @@ subroutine ocn_tracer_hmix_Redi_tend(meshPool, layerThickness, layerThickEdgeMea
real(kind=RKIND), dimension(:, :), intent(in) :: &
layerThickEdgeMean, &!< Input: mean thickness at edge
zMid, &!< Input: Z coordinate at the center of a cell
RediKappaSfcTaper, &!< Input: vertical structure of GM/Redi limiter
RediKappaScaling, &
layerThickness

Expand Down Expand Up @@ -238,14 +237,10 @@ subroutine ocn_tracer_hmix_Redi_tend(meshPool, layerThickness, layerThickEdgeMea

k = minLevelEdgeBot(iEdge)
kappaRediEdgeVal = 0.25_RKIND*(RediKappaScaling(k, cell1) + RediKappaScaling(k, cell2) + &
RediKappaScaling(k + 1, cell1) + RediKappaScaling(k + 1, cell2))* &
0.25_RKIND*(RediKappaSfcTaper(k, cell1) + RediKappaSfcTaper(k, cell2) + &
RediKappaSfcTaper(k + 1, cell1) + RediKappaSfcTaper(k + 1, cell2))
RediKappaScaling(k + 1, cell1) + RediKappaScaling(k + 1, cell2))
k = 1
kappaRediEdgeVal = 0.25_RKIND*(RediKappaScaling(k, cell1) + RediKappaScaling(k, cell2) + &
RediKappaScaling(k + 1, cell1) + RediKappaScaling(k + 1, cell2))* &
0.25_RKIND*(RediKappaSfcTaper(k, cell1) + RediKappaSfcTaper(k, cell2) + &
RediKappaSfcTaper(k + 1, cell1) + RediKappaSfcTaper(k + 1, cell2))
RediKappaScaling(k + 1, cell1) + RediKappaScaling(k + 1, cell2))

do iTr = 1, nTracers
tracer_turb_flux = tracers(iTr, k, cell2) - tracers(iTr, k, cell1)
Expand Down Expand Up @@ -280,9 +275,7 @@ subroutine ocn_tracer_hmix_Redi_tend(meshPool, layerThickness, layerThickEdgeMea
do k = minLevelEdgeBot(iEdge) + 1, maxLevelEdgeTop(iEdge) - 1

kappaRediEdgeVal = 0.25_RKIND*(RediKappaScaling(k, cell1) + RediKappaScaling(k, cell2) + &
RediKappaScaling(k + 1, cell1) + RediKappaScaling(k + 1, cell2))*0.25_RKIND* &
(RediKappaSfcTaper(k, cell1) + RediKappaSfcTaper(k, cell2) + &
RediKappaSfcTaper(k + 1, cell1) + RediKappaSfcTaper(k + 1, cell2))
RediKappaScaling(k + 1, cell1) + RediKappaScaling(k + 1, cell2))

do iTr = 1, nTracers
! \kappa_2 \nabla \phi on edge
Expand Down Expand Up @@ -319,9 +312,7 @@ subroutine ocn_tracer_hmix_Redi_tend(meshPool, layerThickness, layerThickEdgeMea

k = maxLevelEdgeTop(iEdge)
kappaRediEdgeVal = 0.25_RKIND*(RediKappaScaling(k, cell1) + RediKappaScaling(k, cell2) + &
RediKappaScaling(k + 1, cell1) + RediKappaScaling(k + 1, cell2))* &
0.25_RKIND*(RediKappaSfcTaper(k, cell1) + RediKappaSfcTaper(k, cell2) + &
RediKappaSfcTaper(k + 1, cell1) + RediKappaSfcTaper(k + 1, cell2))
RediKappaScaling(k + 1, cell1) + RediKappaScaling(k + 1, cell2))

do iTr = 1, nTracers
tracer_turb_flux = tracers(iTr, k, cell2) - tracers(iTr, k, cell1)
Expand Down Expand Up @@ -376,8 +367,8 @@ subroutine ocn_tracer_hmix_Redi_tend(meshPool, layerThickness, layerThickEdgeMea
! 2.0 in next line is because a dot product on a C-grid
! requires a factor of 1/2 to average to the cell center.
flux_term3 = rediKappaCell(iCell)*2.0_RKIND* &
(RediKappaSfcTaper(k, iCell)*RediKappaScaling(k, iCell)*fluxRediZTop(iTr, k) &
* invAreaCell - RediKappaSfcTaper(k + 1, iCell)*RediKappaScaling(k + 1, iCell) &
(RediKappaScaling(k, iCell)*fluxRediZTop(iTr, k) &
* invAreaCell - RediKappaScaling(k + 1, iCell) &
*fluxRediZTop(iTr, k + 1) * invAreaCell)

redi_term3(iTr, k, iCell) = redi_term3(iTr, k, iCell) + flux_term3
Expand Down Expand Up @@ -446,9 +437,9 @@ subroutine ocn_tracer_hmix_Redi_tend(meshPool, layerThickness, layerThickEdgeMea
do k = minLevelCell(iCell), maxLevelCell(iCell)
do iTr = 1, ntracers
tend(iTr, k, iCell) = tend(iTr, k, iCell) + rediKappaCell(iCell)*2.0_RKIND* &
(RediKappaSfcTaper(k, iCell)*RediKappaScaling(k, iCell)* &
redi_term3_topOfCell(iTr, k, iCell) - RediKappaSfcTaper(k + 1, iCell) &
*RediKappaScaling(k + 1, iCell)*redi_term3_topOfCell(iTr, k + 1, iCell))
(RediKappaScaling(k, iCell)* &
redi_term3_topOfCell(iTr, k, iCell) - &
RediKappaScaling(k + 1, iCell)*redi_term3_topOfCell(iTr, k + 1, iCell))
end do
end do

Expand Down