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

Add flux-corrected transport scheme to MALI #70

Merged
merged 46 commits into from
Sep 6, 2023
+3,956 −44
Merged

Add flux-corrected transport scheme to MALI #70

Show file tree
Hide file tree
Changes from 1 commit
Commits
Show all changes
46 commits
Select commit Hold shift + click to select a range
119cfd9
Add calls to higher order advection scheme
trhille Oct 10, 2022
7532330
Copy higher-order ocean advection routines
trhille Oct 11, 2022
aac672a
Update Makefile to include fct modules
trhille Mar 8, 2023
47d6154
Initialize fct from tracer_setup
trhille Mar 8, 2023
5dd26c6
Add li_mesh module
trhille Mar 9, 2023
e549c55
Comment out irrelevant parameters and routine in li_mesh
trhille Mar 9, 2023
f97dd56
Add li_config module
trhille Mar 9, 2023
ab8ec14
Add fct variables to Registry
trhille Mar 9, 2023
494d092
Remove ocean variables from li_advection_fct and li_advection_fct_shared
trhille Mar 9, 2023
d9f9317
Update shared Makefile to include li_config
trhille Mar 9, 2023
9518f49
Fix call to li_advection_fct_tend
trhille Mar 9, 2023
30d387c
Fix typo in call to li_config
trhille Mar 9, 2023
414eac6
Remove call to tracer_advection_vert_flx
trhille Mar 9, 2023
83f5168
Clean up mesh variable definitions in li_mesh
trhille Mar 10, 2023
79edb95
Fix public parameter allocations
trhille Mar 13, 2023
3c42538
Clean up li_advection_fct_tend
trhille Mar 13, 2023
a364f7a
Add support for fo thickness with fct tracer advection
trhille Mar 14, 2023
e85d3fd
Apply tendencies to tracers with fct advection
trhille Mar 16, 2023
6da1573
Pass layerThicknessOld to fct tracer advection routine
trhille Mar 16, 2023
8734303
Add support for fct thickness advection
trhille Mar 16, 2023
8d08525
Use cellMask_dynamicMargin to calculate advMaskHighOrder
trhille Mar 17, 2023
36ea0a4
Use dynamic ice mask to define advMaskHighOrder
trhille Mar 22, 2023
12981f5
Use boundaryCell to define advMaskHighOrder
trhille Mar 22, 2023
2c66078
Update nTracers as tracers are added
trhille Mar 22, 2023
ab87706
Allocate arrays after nTracers has been calculated
trhille Mar 23, 2023
0f66471
Pass array of 1s as dummy tracer for fct thickness advection
trhille Mar 23, 2023
7e3d897
Fix bug in marking boundaryCell
trhille Mar 23, 2023
0e6fbb9
Add passiveTracer to help verify advection schemes
trhille Mar 23, 2023
697f7f9
Increase max number of tracers to accomodate passiveTracer
trhille Mar 23, 2023
dae96e7
Pass layerThickness as tracer instead of array of 1s
trhille Mar 24, 2023
4c87809
Change normalThicknessFlux, layerThickness, and tracer definitions
trhille Mar 27, 2023
aa85232
Try new mask for 2nd order terms
trhille Apr 19, 2023
70ec7bb
Make 2nd order and 3rd-4th order masks mutually exclusive
trhille Apr 19, 2023
4104681
Clean up 2nd order mask
trhille Apr 20, 2023
4f313f1
Call li_tracer_advection_fct_tend for thickness and tracers separately
trhille Apr 20, 2023
5889ba6
Make fct conserve mass
trhille Apr 21, 2023
fd6792b
Fix first order flux at ice edge
trhille Apr 21, 2023
8d2f023
Make conserve tracer volume
trhille Aug 3, 2023
38f3f2b
Simple cleanup after code review
trhille Aug 18, 2023
7cff4b7
Throw error forr fct thickness with fo tracer
trhille Aug 18, 2023
9530ad5
Make activeTracerHorizontalAdvectionEdgeFlux optional
trhille Aug 18, 2023
aa2999a
Fix bug with activeTracerHorizontalAdvectionEdgeFlux
trhille Aug 20, 2023
35a89ae
Fix the case of fo thickness, none tracer advection
trhille Aug 22, 2023
f43c26b
Only pass layer thickness tracer for first call to fct
trhille Aug 28, 2023
8f218c3
Change passiveTracer to passiveTracer2d
trhille Aug 28, 2023
52d9a0b
Enable fct thickness advection without tracer advection
trhille Sep 1, 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
25 changes: 17 additions & 8 deletions components/mpas-albany-landice/src/mode_forward/mpas_li_advection.F
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -184,8 +184,8 @@ subroutine li_advection_thickness_tracers(&
thermalCellMaskField

! Allocatable arrays need for flux-corrected transport advection
real (kind=RKIND), dimension(:,:), allocatable :: tendLayerThickness
real (kind=RKIND), dimension(:,:,:), allocatable :: tend
real (kind=RKIND), dimension(:,:,:), allocatable :: activeTracerHorizontalAdvectionEdgeFlux

integer, dimension(:), pointer :: &
cellMask, & ! integer bitmask for cells
Expand Down Expand Up @@ -403,10 +403,10 @@ subroutine li_advection_thickness_tracers(&

if ((trim(config_thickness_advection) .eq. 'fct') .or. &
(trim(config_tracer_advection) .eq. 'fct' ) ) then
allocate(tendLayerThickness(nVertLevels,nCells+1))
tendLayerThickness(:,:) = 0.0_RKIND
allocate(tend(nTracers,nVertLevels,nCells+1))
tend(:,:,:) = 0.0_RKIND
allocate(activeTracerHorizontalAdvectionEdgeFlux(nTracers,nVertLevels+1,nEdges+1))
activeTracerHorizontalAdvectionEdgeFlux(:,:,:) = 0.0_RKIND
endif

! Transport thickness and tracers using a first-order upwind scheme
Expand Down Expand Up @@ -469,16 +469,20 @@ subroutine li_advection_thickness_tracers(&
! Reset tracers after fo advection for fct advection
advectedTracers(:,:,:) = advectedTracersOld(:,:,:)

! Pass FO upwind normalThicknessFlux (layerThicknessEdge * layerNormalVelocity)
! to fct tracer routine
call li_tracer_advection_fct_tend(&
tend, advectedTracers, layerThicknessOld, &
layerThicknessEdge * layerNormalVelocity, 0 * normalVelocity, dt, &
nTracers, computeBudgets=.false.)!{{{
nTracers, activeTracerHorizontalAdvectionEdgeFlux, computeBudgets=.false.)!{{{
Copy link

Choose a reason for hiding this comment

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

We may want to move activeTracerHorizontalAdvectionEdgeFlux to Registry so it could be output given its importance in the algorithm. But we'd want to refactor things to avoid it being nTracers big first. So let's wait until we encounter a need to visualize it before worrying about that.

elseif ((trim(config_thickness_advection) .eq. 'fct') .and. &
trim(config_tracer_advection) .eq. 'fct') then
! Call fct routine for thickness first, and use activeTracerHorizontalAdvectionEdgeFlux
! returned by that call as normalThicknessFlux for call to tracer fct
call li_tracer_advection_fct_tend(&
tend, advectedTracers, layerThicknessOld * 0.0_RKIND + 1.0_RKIND, &
matthewhoffman marked this conversation as resolved.
Show resolved Hide resolved
layerNormalVelocity, 0 * normalVelocity, dt, &
nTracers, computeBudgets=.false.)
nTracers, activeTracerHorizontalAdvectionEdgeFlux, computeBudgets=.false.)
! layerThickness is last tracer. However, for some reason
! this: layerThickness(:,:) = advectedTracers(nTracers,:,:) does not conserve mass!
! This does conserve mass:
Copy link

Choose a reason for hiding this comment

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

@dengwirda, do you have any ideas why this might be? We could meet with you to explain the details if no explanation is obvious.

Expand All @@ -487,10 +491,13 @@ subroutine li_advection_thickness_tracers(&
! TODO: determine whether we need to treat other tracer tendencies
Copy link

Choose a reason for hiding this comment

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

@trhille , have you looked into this question?

Copy link
Author

Choose a reason for hiding this comment

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

Yes, it doesn't work for tracers, for some reason.

! the same as layerThickness for conservation.
advectedTracers(:,:,:) = advectedTracersOld(:,:,:)
tend(:,:,:) = 0.0_RKIND
! Call fct for tracers, using activeTracerHorizontalAdvectionEdgeFlux
! from fct thickness advection as normalThicknessFlux
call li_tracer_advection_fct_tend(&
tend, advectedTracers, layerThicknessOld, &
layerThicknessEdge * layerNormalVelocity, 0 * normalVelocity, dt, &
nTracers, computeBudgets=.false.)
activeTracerHorizontalAdvectionEdgeFlux(nTracers,:,:), 0 * normalVelocity, dt, &
trhille marked this conversation as resolved.
Show resolved Hide resolved
nTracers, activeTracerHorizontalAdvectionEdgeFlux, computeBudgets=.false.)
matthewhoffman marked this conversation as resolved.
Show resolved Hide resolved
endif
trhille marked this conversation as resolved.
Show resolved Hide resolved

if (config_print_thickness_advection_info) then
Expand Down Expand Up @@ -690,7 +697,9 @@ subroutine li_advection_thickness_tracers(&
advCoefs, &
advCoefs3rd, &
advMaskHighOrder, &
advMask2ndOrder)
advMask2ndOrder, &
tend, &
activeTracerHorizontalAdvectionEdgeFlux)
endif

! clean up
Expand Down
50 changes: 36 additions & 14 deletions components/mpas-albany-landice/src/mode_forward/mpas_li_advection_fct.F
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -66,15 +66,17 @@ module li_tracer_advection_fct
subroutine li_tracer_advection_fct_tend( &
tend, tracers, layerThickness, &
normalThicknessFlux, w, dt, &
nTracers, computeBudgets)!{{{
nTracers, &
activeTracerHorizontalAdvectionEdgeFlux, &
computeBudgets)!{{{
use li_mesh
!-----------------------------------------------------------------
! Input/Output parameters
!-----------------------------------------------------------------

real (kind=RKIND), dimension(:,:,:), intent(inout) :: &
tend !< [inout] Tracer tendency to which advection added

tend, & !< [inout] Tracer tendency to which advection added
activeTracerHorizontalAdvectionEdgeFlux !< [inout] used to compute higher order normalThicknessFlux
!-----------------------------------------------------------------
! Input parameters
!-----------------------------------------------------------------
Expand Down Expand Up @@ -160,6 +162,21 @@ subroutine li_tracer_advection_fct_tend( &
lowOrderFlx (nVertLevels+1,max(nCells,nEdges)+1), &
highOrderFlx(nVertLevels+1,max(nCells,nEdges)+1))

! Initialize variables so you don't get floating point exceptions
wgtTmp(:) = 0.0_RKIND
flxTmp(:) = 0.0_RKIND
sgnTmp(:) = 0.0_RKIND
tracerCur(:,:) = 0.0_RKIND
tracerMin(:,:) = 0.0_RKIND
tracerMax(:,:) = 0.0_RKIND
hNewInv(:,:) = 0.0_RKIND
hProv(:,:) = 0.0_RKIND
hProvInv(:,:) = 0.0_RKIND
flxIn(:,:) = 0.0_RKIND
flxOut(:,:) = 0.0_RKIND
workTend(:,:) = 0.0_RKIND
lowOrderFlx(:,:) = 0.0_RKIND
highOrderFlx(:,:) = 0.0_RKIND
!$acc enter data create(wgtTmp, flxTmp, sgnTmp, &
!$acc tracerCur, tracerMin, tracerMax, &
!$acc hNewInv, hProv, hProvInv, flxIn, flxOut, &
Expand Down Expand Up @@ -588,17 +605,22 @@ subroutine li_tracer_advection_fct_tend( &
! !$omp parallel
! !$omp do schedule(runtime) private(k)
!#endif
!
! ! TODO: Determine if it is necessary to define activeTracerHorizontalAdvectionEdgeFlux
! !do iEdge = 1,nEdges
! !do k = minLevelEdgeBot(iEdge),maxLevelEdgeTop(iEdge)
! ! ! Save u*h*T flux on edge for analysis. This variable will be
! ! ! divided by h at the end of the time step.
! ! activeTracerHorizontalAdvectionEdgeFlux(iTracer,k,iEdge) = &
! ! (lowOrderFlx(k,iEdge) + highOrderFlx(k,iEdge))/dvEdge(iEdge)
! !enddo
! !enddo
!

! Use activeTracerHorizontalAdvectionEdgeFlux from the call to fct for thickness
! advection as the higher-order normalThicknessFlux for fct tracer advection.
do iEdge = 1,nEdges
do k = 1, nVertLevels
! Save u*h*T flux on edge for analysis. This variable will be
! divided by h at the end of the time step.
if (dvEdge(iEdge) > 0.0_RKIND) then
activeTracerHorizontalAdvectionEdgeFlux(iTracer,k,iEdge) = &
(lowOrderFlx(k,iEdge) + highOrderFlx(k,iEdge))/dvEdge(iEdge)
else
activeTracerHorizontalAdvectionEdgeFlux(iTracer,k,iEdge) = 0.0_RKIND
endif
enddo
enddo

!#ifndef MPAS_OPENACC
! !$omp end do
!#endif
Expand Down