From 399ad2e5d674dd2a7e2d7e171b3443f9a5e3387f Mon Sep 17 00:00:00 2001 From: Jeremy Lilly Date: Thu, 19 Oct 2023 10:28:30 -0700 Subject: [PATCH 01/34] First commit for FB_LTS: Add bare-bones source file and update Makefile --- .../mpas-ocean/src/mode_forward/Makefile | 12 +- .../mode_forward/mpas_ocn_time_integration.F | 10 +- .../mpas_ocn_time_integration_fblts.F | 172 ++++++++++++++++++ 3 files changed, 190 insertions(+), 4 deletions(-) create mode 100644 components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_fblts.F diff --git a/components/mpas-ocean/src/mode_forward/Makefile b/components/mpas-ocean/src/mode_forward/Makefile index 9a96d02fd0da..b3afade5ca52 100644 --- a/components/mpas-ocean/src/mode_forward/Makefile +++ b/components/mpas-ocean/src/mode_forward/Makefile @@ -5,13 +5,18 @@ OBJS = mpas_ocn_forward_mode.o \ mpas_ocn_time_integration_rk4.o \ mpas_ocn_time_integration_si.o \ mpas_ocn_time_integration_split.o \ - mpas_ocn_time_integration_lts.o + mpas_ocn_time_integration_lts.o \ + mpas_ocn_time_integration_fblts.o all: forward_mode forward_mode: $(OBJS) -mpas_ocn_time_integration.o: mpas_ocn_time_integration_rk4.o mpas_ocn_time_integration_si.o mpas_ocn_time_integration_split.o mpas_ocn_time_integration_lts.o +mpas_ocn_time_integration.o: mpas_ocn_time_integration_rk4.o \ + mpas_ocn_time_integration_si.o \ + mpas_ocn_time_integration_split.o \ + mpas_ocn_time_integration_lts.o \ + mpas_ocn_time_integration_fblts.o mpas_ocn_time_integration_rk4.o: @@ -25,7 +30,8 @@ mpas_ocn_forward_mode.o: mpas_ocn_time_integration.o \ mpas_ocn_time_integration_rk4.o \ mpas_ocn_time_integration_si.o \ mpas_ocn_time_integration_split.o \ - mpas_ocn_time_integration_lts.o + mpas_ocn_time_integration_lts.o \ + mpas_ocn_time_integration_fblts.o clean: $(RM) *.o *.mod *.f90 diff --git a/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration.F b/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration.F index 516f12b5157b..138a338a294a 100644 --- a/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration.F +++ b/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration.F @@ -34,6 +34,7 @@ module ocn_time_integration use ocn_time_integration_split use ocn_time_integration_si use ocn_time_integration_lts + use ocn_time_integration_fblts implicit none @@ -70,7 +71,8 @@ module ocn_time_integration timeIntUnsplitExplicit = 2, &! unsplit-explicit timeIntSemiImplicit = 3, &! Semi-implicit timeIntRK4 = 4, &! 4th-order Runge-Kutta - timeIntLTS = 5 ! local time-stepping + timeIntLTS = 5, &! local time-stepping + timeIntFBLTS = 6 ! forward-backward lts !*********************************************************************** @@ -134,6 +136,8 @@ subroutine ocn_timestep(domain, dt, timeStamp)!{{{ call ocn_time_integrator_rk4(domain, dt) case (timeIntLTS) call ocn_time_integrator_lts(domain, dt) + case (timeIntFBLTS) + call ocn_time_integrator_fblts(domain, dt) end select ! compute time since start of simulation, in days @@ -228,6 +232,10 @@ subroutine ocn_timestep_init(domain, dt, err)!{{{ case ('LTS') timeIntegratorChoice = timeIntLTS call ocn_time_integration_lts_init(domain) + + case ('FB_LTS') + timeIntegratorChoice = timeIntFBLTS + call ocn_time_integration_fblts_init(domain) case default diff --git a/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_fblts.F b/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_fblts.F new file mode 100644 index 000000000000..aff1b3733b31 --- /dev/null +++ b/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_fblts.F @@ -0,0 +1,172 @@ +! Copyright (c) 2021, Los Alamos National Security, LLC (LANS) +! and the University Corporation for Atmospheric Research (UCAR). +! +! Unless noted otherwise source code is licensed under the BSD license. +! Additional copyright and license information can be found in the LICENSE file +! distributed with this code, or at http://mpas-dev.github.com/license.html +! +!||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| +! +! ocn_time_integration_lts +! +!> \brief MPAS barotropic ocean LTS Time integration scheme +!> \author Jeremy Lilly +!> \date November 2022 +!> \details +!> This module contains the LTS init routine and the LTS +!> barotropic ocean time integration scheme with splitting +!> on the fast and slow tendency terms. +! +!----------------------------------------------------------------------- + +module ocn_time_integration_fblts + + use mpas_derived_types + use mpas_pool_routines + use mpas_constants + use mpas_dmpar + use mpas_threading + use mpas_vector_reconstruction + use mpas_spline_interpolation + use mpas_timer + + use ocn_constants + use ocn_tendency + use ocn_diagnostics + use ocn_mesh + use ocn_vmix + use ocn_config + use ocn_diagnostics_variables + use ocn_equation_of_state + use ocn_time_average_coupled + use ocn_time_varying_forcing + + implicit none + private + save + + !-------------------------------------------------------------------- + ! + ! Public member functions + ! + !-------------------------------------------------------------------- + public :: ocn_time_integrator_fblts, & + ocn_time_integration_fblts_init + + contains + + +!||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| +! +! ocn_time_integrator_fblts +! +!> \brief MPAS barotropic ocean FB_LTS time integration scheme +!> \author Jeremy Lilly +!> \date October 2023 +!> \details +!> This routine integrates one timestep (dt) using an FB_LTS time +!> integrator with a splitting of the fast and slow tendency terms +! +!----------------------------------------------------------------------- + + subroutine ocn_time_integrator_fblts(domain, dt)!{{{ + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! Advance model state forward in time by the specified time step using + ! a local time stepping scheme with spltting of the fast and slow tendencies + ! + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + implicit none + + !----------------------------------------------------------------- + ! Input variables + !----------------------------------------------------------------- + + real (kind=RKIND), intent(in) :: & + dt !< [in] time step (sec) to move forward + + !----------------------------------------------------------------- + ! Input/output variables + !----------------------------------------------------------------- + + type (domain_type), intent(inout) :: & + domain !< [inout] model state to advance forward + + end subroutine ocn_time_integrator_fblts!}}} + + +!||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| +! +! ocn_time_integration_fblts_init +! +!> \brief Initialization for MPAS ocean FB_LTS time integration scheme +!> \author Giacomo Capodaglio +!> \date October 2023 +!> \details +!> This routine computes the LTS arrays +! +!----------------------------------------------------------------------- + + subroutine ocn_time_integration_fblts_init(domain)!{{{ + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! Sets the LTS arrays + ! + ! Output: LTS arrays are written + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + implicit none + + !----------------------------------------------------------------- + ! Input/output variables + !----------------------------------------------------------------- + + type (domain_type), intent(inout) :: & + domain !< [inout] model state + + end subroutine ocn_time_integration_fblts_init!}}} + + +!||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| +! +! ocn_lts_tends +! +!> \brief MPAS ocean LTS3 Time integration scheme +!> \author Giacomo Capodaglio +!> \date October 2023 +!> \details +!> This routine calculates the tendencies on different LTS regions +! +!----------------------------------------------------------------------- + + subroutine ocn_lts_tends(statePool, LTSPool, tendPool, timeLevelIn, computeOnFineBig, & + computeOnFineSmall, computeOnCoarse, computeOnInterface)!{{{ + + !----------------------------------------------------------------- + ! Input variables + !----------------------------------------------------------------- + + integer, intent(in) :: & + timeLevelIn, & + computeOnFineBig, & + computeOnFineSmall, & + computeOnCoarse, & + computeOnInterface + + type (mpas_pool_type), intent(in) :: & + statePool !< Input: state variables + + type (mpas_pool_type), intent(in) :: & + LTSPool !< Input: LTS data + + !----------------------------------------------------------------- + ! Input/output variables + !----------------------------------------------------------------- + + type (mpas_pool_type), intent(inout) :: & + tendPool !< Input: tendency variables + + end subroutine ocn_lts_tends!}}} + + +end module ocn_time_integration_fblts + From f04af59900cd29cfd14e059c467602221bae1106 Mon Sep 17 00:00:00 2001 From: Jeremy Lilly Date: Thu, 19 Oct 2023 10:31:43 -0700 Subject: [PATCH 02/34] Add fb-weight variables to the registry --- components/mpas-ocean/src/Registry.xml | 16 +++++++++++++++- 1 file changed, 15 insertions(+), 1 deletion(-) diff --git a/components/mpas-ocean/src/Registry.xml b/components/mpas-ocean/src/Registry.xml index 0c454fbca420..64a1b205086b 100644 --- a/components/mpas-ocean/src/Registry.xml +++ b/components/mpas-ocean/src/Registry.xml @@ -202,7 +202,7 @@ /> + + + + + Date: Thu, 19 Oct 2023 11:15:31 -0700 Subject: [PATCH 03/34] Split ocn_lts_tends into seperate subroutines for thickness and velocity --- .../mpas_ocn_time_integration_fblts.F | 56 ++++++++++++++++--- 1 file changed, 49 insertions(+), 7 deletions(-) diff --git a/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_fblts.F b/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_fblts.F index aff1b3733b31..1d92f348b865 100644 --- a/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_fblts.F +++ b/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_fblts.F @@ -128,7 +128,7 @@ end subroutine ocn_time_integration_fblts_init!}}} !||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| ! -! ocn_lts_tends +! ocn_lts_thick_tend ! !> \brief MPAS ocean LTS3 Time integration scheme !> \author Giacomo Capodaglio @@ -138,8 +138,8 @@ end subroutine ocn_time_integration_fblts_init!}}} ! !----------------------------------------------------------------------- - subroutine ocn_lts_tends(statePool, LTSPool, tendPool, timeLevelIn, computeOnFineBig, & - computeOnFineSmall, computeOnCoarse, computeOnInterface)!{{{ + subroutine ocn_lts_thick_tend(statePool, LTSPool, tendPool, timeLevelIn, computeOnFineBig, & + computeOnFineSmall, computeOnCoarse, computeOnInterface)!{{{ !----------------------------------------------------------------- ! Input variables @@ -153,19 +153,61 @@ subroutine ocn_lts_tends(statePool, LTSPool, tendPool, timeLevelIn, computeOnFin computeOnInterface type (mpas_pool_type), intent(in) :: & - statePool !< Input: state variables + statePool !< Input: state variables type (mpas_pool_type), intent(in) :: & - LTSPool !< Input: LTS data + LTSPool !< Input: LTS data !----------------------------------------------------------------- ! Input/output variables !----------------------------------------------------------------- type (mpas_pool_type), intent(inout) :: & - tendPool !< Input: tendency variables + tendPool !< Input/output: tendency variables - end subroutine ocn_lts_tends!}}} + end subroutine ocn_lts_thick_tend!}}} + + +!||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| +! +! ocn_lts_vel_tend +! +!> \brief MPAS ocean LTS3 Time integration scheme +!> \author Giacomo Capodaglio +!> \date October 2023 +!> \details +!> This routine calculates the tendencies on different LTS regions +! +!----------------------------------------------------------------------- + + subroutine ocn_lts_vel_tend(statePool, LTSPool, tendPool, timeLevelIn, computeOnFineBig, & + computeOnFineSmall, computeOnCoarse, computeOnInterface)!{{{ + + !----------------------------------------------------------------- + ! Input variables + !----------------------------------------------------------------- + + integer, intent(in) :: & + timeLevelIn, & + computeOnFineBig, & + computeOnFineSmall, & + computeOnCoarse, & + computeOnInterface + + type (mpas_pool_type), intent(in) :: & + statePool !< Input: state variables + + type (mpas_pool_type), intent(in) :: & + LTSPool !< Input: LTS data + + !----------------------------------------------------------------- + ! Input/output variables + !----------------------------------------------------------------- + + type (mpas_pool_type), intent(inout) :: & + tendPool !< Input/output: tendency variables + + end subroutine ocn_lts_vel_tend!}}} end module ocn_time_integration_fblts From 0a3f96e2823ddaf90f0894f7383471d7f9435ea1 Mon Sep 17 00:00:00 2001 From: Jeremy Lilly Date: Fri, 20 Oct 2023 09:49:14 -0700 Subject: [PATCH 04/34] First pass at seperate lts tends calculations --- .../mpas_ocn_time_integration_fblts.F | 262 +++++++++++++++++- 1 file changed, 259 insertions(+), 3 deletions(-) diff --git a/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_fblts.F b/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_fblts.F index 1d92f348b865..a9bc1ca8bdc1 100644 --- a/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_fblts.F +++ b/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_fblts.F @@ -7,7 +7,7 @@ ! !||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| ! -! ocn_time_integration_lts +! ocn_time_integration_fblts ! !> \brief MPAS barotropic ocean LTS Time integration scheme !> \author Jeremy Lilly @@ -164,7 +164,125 @@ subroutine ocn_lts_thick_tend(statePool, LTSPool, tendPool, timeLevelIn, compute type (mpas_pool_type), intent(inout) :: & tendPool !< Input/output: tendency variables - + + !----------------------------------------------------------------- + ! Local variables + !----------------------------------------------------------------- + + integer, dimension(:,:), pointer :: & + nCellsInLTSRegion + + integer, dimension(:,:,:), pointer :: & + cellsInLTSRegion + + real (kind=RKIND), dimension(:,:), pointer :: & + layerThicknessTend, & + normalVelocity, & + layerThickness + + integer :: & + iEdge, cell1, cell2, k, ie, iRegion, nRegions, & + ic, i, iCell, kmin, kmax + real (kind=RKIND) :: & + invdcEdge, flux + + !----------------------------------------------------------------- + ! Begin routine + !----------------------------------------------------------------- + + nRegions = 2 + + call mpas_pool_get_array(LTSPool, 'cellsInLTSRegion', cellsInLTSRegion) + call mpas_pool_get_array(LTSPool, 'nCellsInLTSRegion', nCellsInLTSRegion) + + call mpas_pool_get_array(tendPool, 'layerThickness', layerThicknessTend) + + call mpas_pool_get_array(statePool, 'normalVelocity', normalVelocity, timeLevelIn) + call mpas_pool_get_array(statePool, 'layerThickness', layerThickness, timeLevelIn) + + ! calculate layerThickEdgeFlux + do iEdge = 1, nEdgesAll + kmin = minLevelEdgeBot(iEdge) + kmax = maxLevelEdgeTop(iEdge) + cell1 = cellsOnEdge(1,iEdge) + cell2 = cellsOnEdge(2,iEdge) + do k = 1,nVertLevels + ! initialize layerThicknessEdgeMean to avoid divide by + ! zero and NaN problems. + layerThickEdgeFlux(k,iEdge) = -1.0e34_RKIND + end do + do k = kmin,kmax + ! central differenced + layerThickEdgeFlux(k,iEdge) = 0.5_RKIND * & + ( layerThickness(k,cell1) & + + layerThickness(k,cell2) ) + end do + end do + + layerThicknessTend(:, :) = 0.0_RKIND + + if (computeOnInterface == 1) then + do iRegion = 1, nRegions + ! thickness tendency + do ic = 1, nCellsInLTSRegion(iRegion,2) + iCell = cellsInLTSRegion(iRegion,2,ic) + do i = 1, nEdgesOnCell(iCell) + iEdge = edgesOnCell(i, iCell) + do k = minLevelEdgeBot(iEdge), maxLevelEdgeTop(iEdge) + flux = normalVelocity(k, iEdge) * dvEdge(iEdge) * layerThickEdgeFlux(k, iEdge) + layerThicknessTend(k, iCell) = layerThicknessTend(k, iCell) & + + edgeSignOnCell(i, iCell) * flux * invAreaCell(iCell) + end do + end do + end do + end do + end if + + if (computeOnFineBig == 1) then + ! thickness tendency + do ic = 1, nCellsInLTSRegion(1,1) + iCell = cellsInLTSRegion(1,1,ic) + do i = 1, nEdgesOnCell(iCell) + iEdge = edgesOnCell(i, iCell) + do k = minLevelEdgeBot(iEdge), maxLevelEdgeTop(iEdge) + flux = normalVelocity(k, iEdge) * dvEdge(iEdge) * layerThickEdgeFlux(k, iEdge) + layerThicknessTend(k, iCell) = layerThicknessTend(k,iCell) & + + edgeSignOnCell(i, iCell) * flux * invAreaCell(iCell) + end do + end do + end do + end if + + if (computeOnFineSmall == 1) then + ! thickness tendency + do ic = 1, nCellsInLTSRegion(1,3) + iCell = cellsInLTSRegion(1,3,ic) + do i = 1, nEdgesOnCell(iCell) + iEdge = edgesOnCell(i, iCell) + do k = minLevelEdgeBot(iEdge), maxLevelEdgeTop(iEdge) + flux = normalVelocity(k, iEdge) * dvEdge(iEdge) * layerThickEdgeFlux(k, iEdge) + layerThicknessTend(k, iCell) = layerThicknessTend(k,iCell) & + + edgeSignOnCell(i, iCell) * flux * invAreaCell(iCell) + end do + end do + end do + end if + + if (computeOnCoarse == 1) then + ! thickness tendency + do ic = 1, nCellsInLTSRegion(2,1) + iCell = cellsInLTSRegion(2,1,ic) + do i = 1, nEdgesOnCell(iCell) + iEdge = edgesOnCell(i, iCell) + do k = minLevelEdgeBot(iEdge), maxLevelEdgeTop(iEdge) + flux = normalVelocity(k, iEdge) * dvEdge(iEdge) * layerThickEdgeFlux(k, iEdge) + layerThicknessTend(k, iCell) = layerThicknessTend(k,iCell) & + + edgeSignOnCell(i, iCell) * flux * invAreaCell(iCell) + end do + end do + end do + end if + end subroutine ocn_lts_thick_tend!}}} @@ -206,7 +324,145 @@ subroutine ocn_lts_vel_tend(statePool, LTSPool, tendPool, timeLevelIn, computeOn type (mpas_pool_type), intent(inout) :: & tendPool !< Input/output: tendency variables - + + !----------------------------------------------------------------- + ! Local variables + !----------------------------------------------------------------- + + integer, dimension(:,:), pointer :: & + nEdgesInLTSRegion + + integer, dimension(:,:,:), pointer :: & + edgesInLTSRegion + + real (kind=RKIND), dimension(:), pointer :: & + ssh + + real (kind=RKIND), dimension(:,:), pointer :: & + normalVelocityTend, & + normalVelocity, & + layerThickness + + integer :: & + iEdge, cell1, cell2, k, ie, iRegion, nRegions, & + ic, i, iCell, kmin, kmax + real (kind=RKIND) :: & + invdcEdge, betaSelfAttrLoad, flux, ssh_sal_on + + !----------------------------------------------------------------- + ! Begin routine + !----------------------------------------------------------------- + + betaSelfAttrLoad = config_self_attraction_and_loading_beta + nRegions = 2 + + call mpas_pool_get_array(LTSPool, 'edgesInLTSRegion', edgesInLTSRegion) + call mpas_pool_get_array(LTSPool, 'nEdgesInLTSRegion', nEdgesInLTSRegion) + + call mpas_pool_get_array(tendPool, 'normalVelocity', normalVelocityTend) + + call mpas_pool_get_array(statePool, 'normalVelocity', normalVelocity, timeLevelIn) + call mpas_pool_get_array(statePool, 'layerThickness', layerThickness, timeLevelIn) + call mpas_pool_get_array(statePool, 'ssh', ssh, timeLevelIn) + + if (config_use_self_attraction_loading) then + ssh_sal_on = 1.0_RKIND + else + ssh_sal_on = 0.0_RKIND + endif + + ! ssh + do iCell = 1, nCellsAll + k = maxLevelCell(iCell) + zTop(k:nVertLevels,iCell) = -bottomDepth(iCell) + layerThickness(k,iCell) + do k = maxLevelCell(iCell)-1, minLevelCell(iCell), -1 + zTop(k,iCell) = zTop(k+1,iCell) + layerThickness(k ,iCell) + end do + ! copy zTop(1,iCell) into sea-surface height array + ssh(iCell) = zTop(minLevelCell(iCell),iCell) + end do + + normalVelocityTend(:,:) = 0.0_RKIND + + ! interface + if (computeOnInterface == 1) then + do iRegion = 1, nRegions + ! velocity tendency + do ie = 1, nEdgesInLTSRegion(iRegion,2) + iEdge = edgesInLTSRegion(iRegion,2,ie) + cell1 = cellsOnEdge(1,iEdge) + cell2 = cellsOnEdge(2,iEdge) + invdcEdge = 1.0_RKIND / dcEdge(iEdge) + kMin = minLevelEdgeBot(iEdge) + kMax = maxLevelEdgeTop(iEdge) + do k=kMin,kMax + normalVelocityTend(k,iEdge) = normalVelocityTend(k,iEdge) & + - edgeMask(k,iEdge) * invdcEdge & + * ( gravity * ( (ssh(cell2) - ssh(cell1)) & + - (1.0_RKIND - ssh_sal_on) * betaSelfAttrLoad & + * (ssh(cell2) - ssh(cell1)) ) ) + end do + end do + end do + end if + + if (computeOnFineBig == 1) then + ! velocity tendency + do ie = 1, nEdgesInLTSRegion(1,1) + iEdge = edgesInLTSRegion(1,1,ie) + cell1 = cellsOnEdge(1,iEdge) + cell2 = cellsOnEdge(2,iEdge) + invdcEdge = 1.0_RKIND / dcEdge(iEdge) + kMin = minLevelEdgeBot(iEdge) + kMax = maxLevelEdgeTop(iEdge) + do k=kMin,kMax + normalVelocityTend(k,iEdge) = normalVelocityTend(k,iEdge) & + - edgeMask(k,iEdge) * invdcEdge & + * ( gravity * ( (ssh(cell2) - ssh(cell1)) & + - (1.0_RKIND - ssh_sal_on) * betaSelfAttrLoad & + * (ssh(cell2) - ssh(cell1)) ) ) + end do + end do + end if + + if (computeOnFineSmall == 1) then + ! velocity tendency + do ie = 1, nEdgesInLTSRegion(1,3) + iEdge = edgesInLTSRegion(1,3,ie) + cell1 = cellsOnEdge(1,iEdge) + cell2 = cellsOnEdge(2,iEdge) + invdcEdge = 1.0_RKIND / dcEdge(iEdge) + kMin = minLevelEdgeBot(iEdge) + kMax = maxLevelEdgeTop(iEdge) + do k=kMin,kMax + normalVelocityTend(k,iEdge) = normalVelocityTend(k,iEdge) & + - edgeMask(k,iEdge) * invdcEdge & + * ( gravity * ( (ssh(cell2) - ssh(cell1)) & + - (1.0_RKIND - ssh_sal_on) * betaSelfAttrLoad & + * (ssh(cell2) - ssh(cell1)) ) ) + end do + end do + end if + + if (computeOnCoarse == 1) then + ! velocity tendency + do ie = 1, nEdgesInLTSRegion(2,1) + iEdge = edgesInLTSRegion(2,1,ie) + cell1 = cellsOnEdge(1,iEdge) + cell2 = cellsOnEdge(2,iEdge) + invdcEdge = 1.0_RKIND / dcEdge(iEdge) + kMin = minLevelEdgeBot(iEdge) + kMax = maxLevelEdgeTop(iEdge) + do k=kMin,kMax + normalVelocityTend(k,iEdge) = normalVelocityTend(k,iEdge) & + - edgeMask(k,iEdge) * invdcEdge & + * ( gravity * ( (ssh(cell2) - ssh(cell1)) & + - (1.0_RKIND - ssh_sal_on) * betaSelfAttrLoad & + * (ssh(cell2) - ssh(cell1)) ) ) + end do + end do + end if + end subroutine ocn_lts_vel_tend!}}} From d7e73c4cf816052aa172e9c8256fa2560a0f00f8 Mon Sep 17 00:00:00 2001 From: Jeremy Lilly Date: Fri, 20 Oct 2023 09:56:58 -0700 Subject: [PATCH 05/34] Update subroutine descriptions --- .../src/mode_forward/mpas_ocn_time_integration_fblts.F | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_fblts.F b/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_fblts.F index a9bc1ca8bdc1..e1eb060f413d 100644 --- a/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_fblts.F +++ b/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_fblts.F @@ -130,7 +130,7 @@ end subroutine ocn_time_integration_fblts_init!}}} ! ! ocn_lts_thick_tend ! -!> \brief MPAS ocean LTS3 Time integration scheme +!> \brief Calculate thickness tendencies for FB_LTS !> \author Giacomo Capodaglio !> \date October 2023 !> \details @@ -290,7 +290,7 @@ end subroutine ocn_lts_thick_tend!}}} ! ! ocn_lts_vel_tend ! -!> \brief MPAS ocean LTS3 Time integration scheme +!> \brief Calculate velocity tendencies for FB_LTS !> \author Giacomo Capodaglio !> \date October 2023 !> \details From cf0d65aa6b91b49c58b5a9797ba31cc65f10321e Mon Sep 17 00:00:00 2001 From: Jeremy Lilly Date: Fri, 20 Oct 2023 12:21:36 -0700 Subject: [PATCH 06/34] Copy init subroutine from LTS3 source to FB_LTS source --- .../mpas_ocn_time_integration_fblts.F | 153 +++++++++++++++++- 1 file changed, 150 insertions(+), 3 deletions(-) diff --git a/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_fblts.F b/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_fblts.F index e1eb060f413d..d1c096dda6b7 100644 --- a/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_fblts.F +++ b/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_fblts.F @@ -108,11 +108,11 @@ end subroutine ocn_time_integrator_fblts!}}} !----------------------------------------------------------------------- subroutine ocn_time_integration_fblts_init(domain)!{{{ - !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! Sets the LTS arrays ! ! Output: LTS arrays are written - !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! implicit none @@ -121,7 +121,154 @@ subroutine ocn_time_integration_fblts_init(domain)!{{{ !----------------------------------------------------------------- type (domain_type), intent(inout) :: & - domain !< [inout] model state + domain !< Input/output: model state + + !----------------------------------------------------------------- + ! Local variables + !----------------------------------------------------------------- + + type (block_type), pointer :: & + block + + type (mpas_pool_type), pointer :: & + LTSPool, & + meshPool + + integer, dimension(:), allocatable :: & + isLTSRegionEdgeAssigned + + integer, dimension(:), pointer :: & + LTSRegion + + integer, dimension(:,:), pointer :: & + nCellsInLTSRegion, & + nEdgesInLTSRegion + + integer, dimension(:,:,:), pointer :: & + cellsInLTSRegion, & + edgesInLTSRegion + + integer, dimension(2) :: & + minMaxLTSRegion + + integer :: & + i, iCell, iEdge, iRegion, coarseRegions, & + fineRegions, fineRegionsM1 + + !----------------------------------------------------------------- + ! Begin routine + !----------------------------------------------------------------- + + minMaxLTSRegion(1) = 1 + minMaxLTSRegion(2) = 2 + + block => domain % blocklist + call mpas_pool_get_subpool(block % structs, 'LTS', LTSPool) + call mpas_pool_get_subpool(block % structs, ',mesh', meshPool) + + call mpas_pool_get_array(LTSPool, 'LTSRegion', LTSRegion) + call mpas_pool_get_array(LTSPool, 'cellsInLTSRegion', cellsInLTSRegion) + call mpas_pool_get_array(LTSPool, 'nCellsInLTSRegion', nCellsInLTSRegion) + call mpas_pool_get_array(LTSPool, 'edgesInLTSRegion', edgesInLTSRegion) + call mpas_pool_get_array(LTSPool, 'nEdgesInLTSRegion', nEdgesInLTSRegion) + + ! LTS Regions code: + ! 1 = fine + ! 2 = coarse + ! 3 = interface layer 1 + ! 4 = interface layer 2 + ! 5 = fine (to advance when doing 1st, 2nd, 3rd stage on interface) + + nCellsInLTSRegion(:,:) = 0 + nEdgesInLTSRegion(:,:) = 0 + + ! This is a loop to build the lists of elements in the fine, coarse, + ! and interface regions. Only loops up to nCellsOwned because in the + ! time-stepping we only want to advance the cells owned by the MPI process + do iCell = 1, nCellsOwned + do iRegion = 1,2 + if (iRegion == minMaxLTSRegion(iRegion)) then + if(LTSRegion(iCell) == minMaxLTSRegion(iRegion)) then + nCellsInLTSRegion(iRegion,1) = nCellsInLTSRegion(iRegion,1) + 1 + cellsInLTSRegion(iRegion,1,nCellsInLTSRegion(iRegion,1)) = iCell + end if + if(LTSRegion(iCell) == (minMaxLTSRegion(iRegion) + 2) ) then + nCellsInLTSRegion(iRegion,2) = nCellsInLTSRegion(iRegion,2) + 1 + cellsInLTSRegion(iRegion,2,nCellsInLTSRegion(iRegion,2)) = iCell + end if + end if + end do + if (LTSRegion(iCell) == 5) then + nCellsInLTSRegion(1,3) = nCellsInLTSRegion(1,3) + 1 + cellsInLTSRegion(1,3,nCellsInLTSRegion(1,3)) = iCell + end if + end do + + ! Below we fill out the lists for the edges, according to the LTSRegion + ! that have been assigned to the cells. We move from the fine to the + ! coarse (i.e. from the fine to the nearest LTS region in the direction + ! of the coarse). Note that edges shared between cells of different LTS + ! regions are owned by the cell in the LTS region closest to the fine + ! region, see Figure 3 in "Conservative explicit local time-stepping + ! schemes for the shallow water equations" by Hoang et al. (halo edges + ! however are owned by whatever processor they are initially assigned to). + + allocate(isLTSRegionEdgeAssigned(nEdgesOwned)) + isLTSRegionEdgeAssigned(:) = 0 + + do iCell = 1, nCellsInLTSRegion(1,1) + do i = 1, nEdgesOnCell(cellsInLTSRegion(1,1,iCell)) + iEdge = edgesOnCell(i,cellsInLTSRegion(1,1,iCell)) + if (iEdge .le. nEdgesOwned) then + if (isLTSRegionEdgeAssigned(iEdge) == 0) then + nEdgesInLTSRegion(1,1) = nEdgesInLTSRegion(1,1) + 1 + edgesInLTSRegion(1,1, nEdgesInLTSRegion(1,1)) = iEdge + isLTSRegionEdgeAssigned(iEdge) = 1 + end if + end if + end do + end do + + fineRegions = 3 + fineRegionsM1 = 2 + do iRegion = 1, fineRegionsM1 + do iCell = 1, nCellsInLTSRegion(1, fineRegions - iRegion + 1) + do i = 1, nEdgesOnCell(cellsInLTSRegion(1, fineRegions - iRegion + 1, iCell)) + iEdge = edgesOnCell(i,cellsInLTSRegion(1, fineRegions - iRegion + 1, iCell)) + if (iEdge .le. nEdgesOwned) then + if (isLTSRegionEdgeAssigned(iEdge) == 0) then + nEdgesInLTSRegion(1, fineRegions - iRegion + 1) & + = nEdgesInLTSRegion(1, fineRegions - iRegion + 1) + 1 + edgesInLTSRegion(1, fineRegions - iRegion + 1, & + nEdgesInLTSRegion(1, fineRegions - iRegion + 1)) & + = iEdge + isLTSRegionEdgeAssigned(iEdge) = 1 + end if + end if + end do + end do + end do + + coarseRegions = 2 + do iRegion = 1, coarseRegions + do iCell = 1, nCellsInLTSRegion(2, coarseRegions - iRegion + 1) + do i = 1, nEdgesOnCell(cellsInLTSRegion(2,coarseRegions - iRegion + 1,iCell)) + iEdge = edgesOnCell(i,cellsInLTSRegion(2,coarseRegions - iRegion + 1,iCell)) + if (iEdge .le. nEdgesOwned) then + if (isLTSRegionEdgeAssigned(iEdge) == 0) then + nEdgesInLTSRegion(2, coarseRegions - iRegion + 1) & + = nEdgesInLTSRegion(2, coarseRegions - iRegion + 1) + 1 + edgesInLTSRegion(2, coarseRegions - iRegion + 1, & + nEdgesInLTSRegion(2, coarseRegions - iRegion + 1)) & + = iEdge + isLTSRegionEdgeAssigned(iEdge) = 1 + end if + end if + end do + end do + end do + + deallocate(isLTSRegionEdgeAssigned) end subroutine ocn_time_integration_fblts_init!}}} From 231e5817b7a142acef827eac49b50b25e54645b3 Mon Sep 17 00:00:00 2001 From: Jeremy Lilly Date: Mon, 23 Oct 2023 10:07:47 -0700 Subject: [PATCH 07/34] Remove unneeded mesh pool from init --- .../src/mode_forward/mpas_ocn_time_integration_fblts.F | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_fblts.F b/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_fblts.F index d1c096dda6b7..6d8bfe23d78b 100644 --- a/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_fblts.F +++ b/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_fblts.F @@ -131,8 +131,7 @@ subroutine ocn_time_integration_fblts_init(domain)!{{{ block type (mpas_pool_type), pointer :: & - LTSPool, & - meshPool + LTSPool integer, dimension(:), allocatable :: & isLTSRegionEdgeAssigned @@ -164,7 +163,6 @@ subroutine ocn_time_integration_fblts_init(domain)!{{{ block => domain % blocklist call mpas_pool_get_subpool(block % structs, 'LTS', LTSPool) - call mpas_pool_get_subpool(block % structs, ',mesh', meshPool) call mpas_pool_get_array(LTSPool, 'LTSRegion', LTSRegion) call mpas_pool_get_array(LTSPool, 'cellsInLTSRegion', cellsInLTSRegion) From 674b08bae57e1f4c363a09ef2b9e7dc8e3dbab01 Mon Sep 17 00:00:00 2001 From: Jeremy Lilly Date: Mon, 23 Oct 2023 11:31:07 -0700 Subject: [PATCH 08/34] Add lts prep step --- .../mpas_ocn_time_integration_fblts.F | 288 +++++++++++++++++- 1 file changed, 284 insertions(+), 4 deletions(-) diff --git a/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_fblts.F b/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_fblts.F index 6d8bfe23d78b..bc03de7b265b 100644 --- a/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_fblts.F +++ b/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_fblts.F @@ -70,11 +70,12 @@ module ocn_time_integration_fblts !----------------------------------------------------------------------- subroutine ocn_time_integrator_fblts(domain, dt)!{{{ - !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - ! Advance model state forward in time by the specified time step using - ! a local time stepping scheme with spltting of the fast and slow tendencies + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! Advance model state forward in time by the specified time step + ! using a local time stepping scheme with spltting of the fast and + ! slow tendencies ! - !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! implicit none @@ -92,6 +93,285 @@ subroutine ocn_time_integrator_fblts(domain, dt)!{{{ type (domain_type), intent(inout) :: & domain !< [inout] model state to advance forward + !----------------------------------------------------------------- + ! Local variables + !----------------------------------------------------------------- + + integer :: & + iCell, iEdge, iRegion, k, ic, ie, im, & ! iterators + M, & ! M = dtCoarse / dtFine + nRegions ! number of interface regions (two) + + type (block_type), pointer :: & + block ! structure with subdomain data + + type (mpas_pool_type), pointer :: & + tendPool, & ! structure holding tendencies + statePool, & ! structure holding state variables + meshPool, & ! structure holding mesh variables + verticalMeshPool, & ! structure holding mesh variables + forcingPool, & ! structure holding forcing variables + scratchPool, & ! structure holding temporary variables + tracersPool ! structure holding tracers variables + + ! LTS Pools + type (mpas_pool_type), pointer :: & + LTSPool, & ! structure holding LTS variables + tendSlowPool, & ! structure holding the slow tendency variables + tendSum1stPool, & ! structure holding one of the correction terms for the interface + tendSum2ndPool, & ! structure holding one of the correction terms for the interface + tendSum3rdPool, & ! structure holding one of the correction terms for the interface + prevTendSlowPool, nextTendSlowPool, & ! structures containing intermediate data + prevTendSum1stPool, nextTendSum1stPool, & ! structures containing intermediate data + prevTendSum2ndPool, nextTendSum2ndPool, & ! structures containing intermediate data + prevTendSum3rdPool, nextTendSum3rdPool ! structures containing intermediate data + + ! Tend Array Pointers + real (kind=RKIND), dimension(:,:), pointer :: & + normalVelocityTend, & ! normal velocity fast tendency + layerThicknessTend, & ! layer thickness tendency + normalVelocityTendSlow, & ! normal velocity slow tendency + normalVelocityTendSum1st, & ! one of the normal velocity correction terms for the interface + layerThicknessTendSum1st, & ! one of the layer thickness correction terms for the interface + normalVelocityTendSum2nd, & ! one of the normal velocity correction terms for the interface + layerThicknessTendSum2nd, & ! one of the layer thickness correction terms for the interface + normalVelocityTendSum3rd, & ! one of the normal velocity correction terms for the interface + layerThicknessTendSum3rd ! one of the layer thickness correction terms for the interface + + ! State Array Pointers + real (kind=RKIND), dimension(:,:), pointer :: & + normalVelocityCur, & ! normal velocity at time n + normalVelocityNew, & ! normal velocity at time n+1 + normalVelocityFirstStage, & ! normal velocity at first stage of LTS + normalVelocitySecondStage, & ! normal velocity at second stage of LTS + layerThicknessCur, & ! layer thickness at time n + layerThicknessNew, & ! layer thickness at time n+1 + layerThicknessFirstStage, & ! layer thickness at first stage of LTS + layerThicknessSecondStage ! layer thickness at second stage of LTS + + ! LTS objects + real (kind=RKIND) :: & + dtFine, & ! fine dt, defined as dt / M + alpha, alphaHat, beta, betaHat, gam, gamHat ! interpolation coefficients for the predictor step of LTS + integer, dimension(:,:), pointer :: & + nCellsInLTSRegion, & ! number of cells in a given LTS region + nEdgesInLTSRegion ! number of edges in a given LTS region + integer, dimension(:,:,:), pointer :: & + cellsInLTSRegion, & ! list of cells in a given LTS region + edgesInLTSRegion ! list of edges in a given LTS region + real (kind=RKIND) :: & + weight1st, weight2nd, weight3rd ! coefficents each RK stage + real (kind=RKIND) :: & + weightTendSum1st, weightTendSum2nd, weightTendSum3rd ! coefficients for the interface correction + + integer err + err = 0 + + !----------------------------------------------------------------- + ! Begin routine + !----------------------------------------------------------------- + + if (.not. config_disable_tr_all_tend) then + call mpas_log_write("ERROR: tracers are not currently implemented for local time-stepping") + call mpas_log_write("config_disable_tr_all_tend should be true in the namelist file") + call abort + end if + + + call mpas_timer_start("FB_LTS time-step prep") + + ! LTS parameters + M = config_dt_scaling_LTS + nRegions = 2 + dtFine = dt / M + + ! Weights for RK stages + weight1st = 1.0_RKIND / 3.0_RKIND + weight2nd = 1.0_RKIND / 2.0_RKIND + weight3rd = 1.0_RKIND + + block => domain % blocklist + + ! Retrieve model state, pools + call mpas_pool_get_subpool(block%structs, 'mesh', meshPool) + call mpas_pool_get_subpool(block%structs, 'verticalMesh', verticalMeshPool) + call mpas_pool_get_subpool(block%structs, 'state', statePool) + call mpas_pool_get_subpool(block%structs, 'forcing', forcingPool) + call mpas_pool_get_subpool(block%structs, 'LTS', LTSPool) + call mpas_pool_get_subpool(block%structs, 'tend', tendPool) + call mpas_pool_get_subpool(block%structs, 'scratch', scratchPool) + call mpas_pool_get_subpool(statePool, 'tracers', tracersPool) + + ! Retrieve state variables at necessary time levels + call mpas_pool_get_array(statePool, 'normalVelocity', & + normalVelocityCur, 1) + call mpas_pool_get_array(statePool, 'normalVelocity', & + normalVelocityNew, 2) + call mpas_pool_get_array(statePool, 'normalVelocity', & + normalVelocityFirstStage, 3) + call mpas_pool_get_array(statePool, 'normalVelocity', & + normalVelocitySecondStage, 4) + call mpas_pool_get_array(statePool, 'layerThickness', & + layerThicknessCur, 1) + call mpas_pool_get_array(statePool, 'layerThickness', & + layerThicknessNew, 2) + call mpas_pool_get_array(statePool, 'layerThickness', & + layerThicknessFirstStage, 3) + call mpas_pool_get_array(statePool, 'layerThickness', & + layerThicknessSecondStage, 4) + + ! Retrieve tendency variables + call mpas_pool_get_array(tendPool, 'normalVelocity', normalVelocityTend) + call mpas_pool_get_array(tendPool, 'layerThickness', layerThicknessTend) + + ! Retrieve LTS arrays + call mpas_pool_get_array(LTSPool, 'cellsInLTSRegion', cellsInLTSRegion) + call mpas_pool_get_array(LTSPool, 'nCellsInLTSRegion', nCellsInLTSRegion) + call mpas_pool_get_array(LTSPool, 'edgesInLTSRegion', edgesInLTSRegion) + call mpas_pool_get_array(LTSPool, 'nEdgesInLTSRegion', nEdgesInLTSRegion) + + ! Create and retrieve additional pools for LTS + !call mpas_pool_create_pool(tendSum1stPool) + !call mpas_pool_clone_pool(tendPool, tendSum1stPool, 1) + !call mpas_pool_create_pool(tendSum2ndPool) + !call mpas_pool_clone_pool(tendPool, tendSum2ndPool, 1) + call mpas_pool_create_pool(tendSum3rdPool) + call mpas_pool_clone_pool(tendPool, tendSum3rdPool, 1) + call mpas_pool_create_pool(tendSlowPool) + call mpas_pool_clone_pool(tendPool, tendSlowPool, 1) + + !call mpas_pool_add_subpool(block % structs, 'tend_sum_1st', tendSum1stPool) + !call mpas_pool_add_subpool(block % structs, 'tend_sum_2nd', tendSum2ndPool) + call mpas_pool_add_subpool(block % structs, 'tend_sum_3rd', tendSum3rdPool) + call mpas_pool_add_subpool(block % structs, 'tend_slow', tendSlowPool) + + call mpas_pool_get_array(tendSlowPool, 'normalVelocity', & + normalVelocityTendSlow) + !call mpas_pool_get_array(tendSum1stPool, 'normalVelocity', & + ! normalVelocityTendSum1st) + !call mpas_pool_get_array(tendSum1stPool, 'layerThickness', & + ! layerThicknessTendSum1st) + !call mpas_pool_get_array(tendSum2ndPool, 'normalVelocity', & + ! normalVelocityTendSum2nd) + !call mpas_pool_get_array(tendSum2ndPool, 'layerThickness', & + ! layerThicknessTendSum2nd) + call mpas_pool_get_array(tendSum3rdPool, 'normalVelocity', & + normalVelocityTendSum3rd) + call mpas_pool_get_array(tendSum3rdPool, 'layerThickness', & + layerThicknessTendSum3rd) + + ! Init variables + do iEdge = 1, nEdgesAll ! do we need this? + do k = 1, maxLevelEdgeTop(iEdge) + normalVelocityNew(k, iEdge) = normalVelocityCur(k, iEdge) + normalVelocityFirstStage(k, iEdge) = normalVelocityCur(k, iEdge) + normalVelocitySecondStage(k, iEdge) = normalVelocityCur(k, iEdge) + end do + end do + + do iCell = 1, nCellsAll ! do we need this? + do k = 1, maxLevelCell(iCell) + layerThicknessNew(k, iCell) = layerThicknessCur(k, iCell) + layerThicknessFirstStage(k, iCell) = layerThicknessCur(k, iCell) + layerThicknessSecondStage(k, iCell) = layerThicknessCur(k, iCell) + end do + end do + + !normalVelocityTendSum1st(:,:) = 0.0_RKIND + !layerThicknessTendSum1st(:,:) = 0.0_RKIND + + !normalVelocityTendSum2nd(:,:) = 0.0_RKIND + !layerThicknessTendSum2nd(:,:) = 0.0_RKIND + + normalVelocityTendSum3rd(:,:) = 0.0_RKIND + layerThicknessTendSum3rd(:,:) = 0.0_RKIND + + if (associated(block % prev)) then + !call mpas_pool_get_subpool(block % prev % structs, 'tend_sum_1st', tendSum1stPool) + !call mpas_pool_get_subpool(block % prev % structs, 'tend_sum_2nd', tendSum2ndPool) + call mpas_pool_get_subpool(block % prev % structs, 'tend_sum_3rd', tendSum3rdPool) + call mpas_pool_get_subpool(block % prev % structs, 'tend_slow', tendSlowPool) + else + !nullify(prevTendSum1stPool) + !nullify(prevTendSum2ndPool) + nullify(prevTendSum3rdPool) + nullify(prevTendSlowPool) + end if + + if (associated(block % next)) then + !call mpas_pool_get_subpool(block % next % structs, 'tend_sum_1st', nextTendSum1stPool) + !call mpas_pool_get_subpool(block % next % structs, 'tend_sum_2nd', nextTendSum2ndPool) + call mpas_pool_get_subpool(block % next % structs, 'tend_sum_3rd', nextTendSum3rdPool) + call mpas_pool_get_subpool(block % next % structs, 'tend_slow', nextTendSlowPool) + else + !nullify(nextTendSum1stPool) + !nullify(nextTendSum2ndPool) + nullify(nextTendSum3rdPool) + nullify(nextTendSlowPool) + end if + + !call mpas_pool_get_subpool(block % structs, 'tend_sum_1st', tendSum1stPool) + !call mpas_pool_get_subpool(block % structs, 'tend_sum_2nd', tendSum2ndPool) + call mpas_pool_get_subpool(block % structs, 'tend_sum_3rd', tendSum3rdPool) + call mpas_pool_get_subpool(block % structs, 'tend_slow', tendSlowPool) + + !if (associated(prevTendSum1stPool) .and. associated(nextTendSum1stPool)) then + ! call mpas_pool_link_pools(tendSum1stPool, prevTendSum1stPool, nextTendSum1stPool) + !else if (associated(prevTendSum1stPool)) then + ! call mpas_pool_link_pools(tendSum1stPool, prevTendSum1stPool) + !else if (associated(nextTendSum1stPool)) then + ! call mpas_pool_link_pools(tendSum1stPool,nextPool=nextTendSum1stPool) + !else + ! call mpas_pool_link_pools(tendSum1stPool) + !end if + + !if (associated(prevTendSum2ndPool) .and. associated(nextTendSum2ndPool)) then + ! call mpas_pool_link_pools(tendSum2ndPool, prevTendSum2ndPool, nextTendSum2ndPool) + !else if (associated(prevTendSum2ndPool)) then + ! call mpas_pool_link_pools(tendSum2ndPool, prevTendSum2ndPool) + !else if (associated(nextTendSum2ndPool)) then + ! call mpas_pool_link_pools(tendSum2ndPool,nextPool=nextTendSum2ndPool) + !else + ! call mpas_pool_link_pools(tendSum2ndPool) + !end if + + if (associated(prevTendSum3rdPool) .and. associated(nextTendSum3rdPool)) then + call mpas_pool_link_pools(tendSum3rdPool, prevTendSum3rdPool, nextTendSum3rdPool) + else if (associated(prevTendSum3rdPool)) then + call mpas_pool_link_pools(tendSum3rdPool, prevTendSum3rdPool) + else if (associated(nextTendSum3rdPool)) then + call mpas_pool_link_pools(tendSum3rdPool,nextPool=nextTendSum3rdPool) + else + call mpas_pool_link_pools(tendSum3rdPool) + end if + + if (associated(prevTendSlowPool) .and. associated(nextTendSlowPool)) then + call mpas_pool_link_pools(tendSlowPool, prevTendSlowPool, nextTendSlowPool) + else if (associated(prevTendSlowPool)) then + call mpas_pool_link_pools(tendSlowPool, prevTendSlowPool) + else if (associated(nextTendSlowPool)) then + call mpas_pool_link_pools(tendSlowPool,nextPool=nextTendSlowPool) + else + call mpas_pool_link_pools(tendSlowPool) + end if + + call mpas_pool_link_parinfo(block, tendSum1stPool) + call mpas_pool_link_parinfo(block, tendSum2ndPool) + call mpas_pool_link_parinfo(block, tendSum3rdPool) + call mpas_pool_link_parinfo(block, tendSlowPool) + + call mpas_timer_stop("FB_LTS time-step prep") + + + call mpas_timer_start("FB_LTS main loop") + + call mpas_timer_stop("FB_LTS main loop") + + + call mpas_timer_start("FB_LTS cleanup phase") + + call mpas_timer_stop("FB_LTS cleanup phase") + end subroutine ocn_time_integrator_fblts!}}} From e00be9089c63cf547ee8d04c24d2416c1fcf12da Mon Sep 17 00:00:00 2001 From: Jeremy Lilly Date: Mon, 23 Oct 2023 15:08:34 -0700 Subject: [PATCH 09/34] Add comments outlining the main steps for FB_LTS --- .../mpas_ocn_time_integration_fblts.F | 161 +++++++++++++++++- 1 file changed, 159 insertions(+), 2 deletions(-) diff --git a/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_fblts.F b/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_fblts.F index bc03de7b265b..9456f7f0f9a3 100644 --- a/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_fblts.F +++ b/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_fblts.F @@ -355,18 +355,175 @@ subroutine ocn_time_integrator_fblts(domain, dt)!{{{ call mpas_pool_link_pools(tendSlowPool) end if - call mpas_pool_link_parinfo(block, tendSum1stPool) - call mpas_pool_link_parinfo(block, tendSum2ndPool) + !call mpas_pool_link_parinfo(block, tendSum1stPool) + !call mpas_pool_link_parinfo(block, tendSum2ndPool) call mpas_pool_link_parinfo(block, tendSum3rdPool) call mpas_pool_link_parinfo(block, tendSlowPool) call mpas_timer_stop("FB_LTS time-step prep") + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! + ! BEGIN FB_LTS SCHEME + ! + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + call mpas_timer_start("FB_LTS main loop") + + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! BEGIN: Slow tendency calculation + ! Calculate the slow tendencies on all LTS regions (only the + ! momentum equation contains slow terms) + call mpas_timer_start("FB_LTS slow tend computation") + + call ocn_tend_vel(domain, tendSlowPool, statePool, forcingPool, 1, & + domain % dminfo, dt) + + call mpas_timer_stop("FB_LTS slow tend computation") + ! END: Slow tendency calculation + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! BEGIN COARSE ADVANCEMENT + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! BEGIN: Thickness stage 1 + ! Compute the first stage of the thickness on fine*, interface 1, + ! interface 2, and coarse + + ! END: Thickness stage 1 + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! BEGIN: Normal velocity stage 1 + ! Compute the first stage of the normal velocity on fine*, + ! interface 1, interface 2, and coarse + + ! END: Normal velocity stage 1 + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! BEGIN: Thickness stage 2 + ! Compute the second stage of the thickness on fine*, interface 1, + ! interface 2, and coarse + + ! END: Thickness stage 2 + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! BEGIN: Normal velocity stage 2 + ! Compute the second stage of the normal velocity on fine*, + ! interface 1, interface 2, and coarse + + ! END: Normal velocity stage 2 + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! BEGIN: Thickness stage 3 + ! Compute the third stage of the thickness on fine*, interface 1, + ! interface 2, and coarse + + ! END: Thickness stage 3 + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! BEGIN: Normal velocity stage 3 + ! Compute the third stage of the normal velocity on, interface 1, + ! and coarse + + ! END: Normal velocity stage 3 + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! END COARSE ADVANCEMENT + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! BEGIN FINE ADVANCEMENT + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + do im = 1,M + + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! BEGIN: Thickness stage 1 + ! Compute the first stage of the thickness on fine + + ! END: Thickness stage 1 + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! BEGIN: Normal velocity stage 1 + ! Compute the first stage of the normal velocity on fine + + ! END: Normal velocity stage 1 + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! BEGIN: Thickness stage 2 + ! Compute the second stage of the thickness on fine + + ! END: Thickness stage 2 + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! BEGIN: Normal velocity stage 2 + ! Compute the second stage of the normal velocity on fine + + ! END: Normal velocity stage 2 + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! BEGIN: Thickness stage 3 + ! Compute the third stage of the thickness on fine and + ! increment the third stage interface correction terms + + ! END: Thickness stage 3 + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! BEGIN: Normal velocity stage 3 + ! Compute the third stage of the normal velocity on fine and + ! increment the thrid stage interface correction terms + + ! END: Normal velocity stage 3 + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + end do + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! END FINE ADVANCEMENT + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! BEGIN INTERFACE CORRECTION + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! END INTERFACE CORRECTION + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! call mpas_timer_stop("FB_LTS main loop") + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! + ! END FB_LTS SCHEME + ! + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + call mpas_timer_start("FB_LTS cleanup phase") From 3d3db37e2cee795a89420476a215e62dd272ee2d Mon Sep 17 00:00:00 2001 From: Jeremy Lilly Date: Mon, 23 Oct 2023 16:47:14 -0700 Subject: [PATCH 10/34] Add first RK stage in coarse cells, the FB average of thickness data for velocity advancement needs to be implemented --- .../mpas_ocn_time_integration_fblts.F | 109 +++++++++++++++++- 1 file changed, 106 insertions(+), 3 deletions(-) diff --git a/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_fblts.F b/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_fblts.F index 9456f7f0f9a3..5d6ccb484dc9 100644 --- a/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_fblts.F +++ b/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_fblts.F @@ -375,12 +375,12 @@ subroutine ocn_time_integrator_fblts(domain, dt)!{{{ ! BEGIN: Slow tendency calculation ! Calculate the slow tendencies on all LTS regions (only the ! momentum equation contains slow terms) - call mpas_timer_start("FB_LTS slow tend computation") + call mpas_timer_start("FB_LTS compute slow tendencies") call ocn_tend_vel(domain, tendSlowPool, statePool, forcingPool, 1, & domain % dminfo, dt) - call mpas_timer_stop("FB_LTS slow tend computation") + call mpas_timer_stop("FB_LTS compute slow tendencies") ! END: Slow tendency calculation !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! @@ -389,11 +389,65 @@ subroutine ocn_time_integrator_fblts(domain, dt)!{{{ ! BEGIN COARSE ADVANCEMENT !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! Update halos for diagnostic variables + call mpas_timer_start("FB_LTS diagnostic halo update") + + call mpas_dmpar_field_halo_exch(domain, 'normalizedRelativeVorticityEdge') + if (config_mom_del4 > 0.0_RKIND) then + call mpas_dmpar_field_halo_exch(domain, 'divergence') + call mpas_dmpar_field_halo_exch(domain, 'relativeVorticity') + end if + + call mpas_timer_stop("FB_LTS diagnostic halo update") + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! BEGIN: Thickness stage 1 ! Compute the first stage of the thickness on fine*, interface 1, ! interface 2, and coarse - + + ! Compute fast tendencies for thickness + call mpas_timer_start("FB_LTS compute fast tendencies") + call mpas_timer_start("FB_LTS compute fast thickness tendencies") + + call ocn_lts_thick_tend(statePool, LTSPool, tendPool, 1, 0, 1, 1, 1) + + call mpas_timer_stop("FB_LTS compute fast thickness tendencies") + call mpas_timer_stop("FB_LTS compute fast tendencies") + + ! Advance thickness solution with first stage + call mpas_timer_start("FB_LTS advance solution") + call mpas_timer_start("FB_LTS advance thickness solution") + + ! Interface layers + do iRegion =1,nRegions + do ic = 1, nCellsInLTSRegion(iRegion,2) + iCell = cellsInLTSRegion(iRegion,2,ic) + do k = 1, maxLevelCell(iCell) + layerThicknessFirstStage(k,iCell) = layerThicknessCur(k,iCell) & + + weight1st * dt * layerThicknessTend(k,iCell) + end do + end do + end do + ! Coarse + do ic = 1, nCellsInLTSRegion(2,1) + iCell = cellsInLTSRegion(2,1,ic) + do k = 1, maxLevelCell(iCell) + layerThicknessFirstStage(k,iCell) = layerThicknessCur(k,iCell) & + + weight1st * dt * layerThicknessTend(k,iCell) + end do + end do + ! Fine layers close to interface layers + do ic = 1, nCellsInLTSRegion(1,3) + iCell = cellsInLTSRegion(1,3,ic) + do k = 1, maxLevelCell(iCell) + layerThicknessFirstStage(k,iCell) = layerThicknessCur(k,iCell) & + + weight1st * dt * layerThicknessTend(k,iCell) + end do + end do + + call mpas_timer_stop("FB_LTS advance thickness solution") + call mpas_timer_stop("FB_LTS advance solution") + ! END: Thickness stage 1 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! @@ -403,9 +457,58 @@ subroutine ocn_time_integrator_fblts(domain, dt)!{{{ ! Compute the first stage of the normal velocity on fine*, ! interface 1, interface 2, and coarse + ! Compute fast tendencies for thickness + call mpas_timer_start("FB_LTS compute fast tendencies") + call mpas_timer_start("FB_LTS compute fast velocity tendencies") + + ! NEED TO CALCULATE h^* + call ocn_lts_vel_tend(statePool, LTSPool, tendPool, 1, 0, 1, 1, 1) + + call mpas_timer_stop("FB_LTS compute fast velocity tendencies") + call mpas_timer_stop("FB_LTS compute fast tendencies") + + ! Advance normal velocity solution with first stage + call mpas_timer_start("FB_LTS advance solution") + call mpas_timer_start("FB_LTS advance velocity solution") + + normalVelocityTend(:,:) = normalVelocityTend(:,:) & + + normalVelocityTendSlow(:,:) + + ! Interface regions + do iRegion =1,nRegions + do ie = 1, nEdgesInLTSRegion(iRegion,2) + iEdge = edgesInLTSRegion(iRegion,2,ie) + normalVelocityFirstStage(:,iEdge) = normalVelocityCur(:,iEdge) & + + weight1st * dt * normalVelocityTend(:,iEdge) + end do + end do + ! Coarse + do ie = 1, nEdgesInLTSRegion(2,1) + iEdge = edgesInLTSRegion(2,1,ie) + normalVelocityFirstStage(:,iEdge) = normalVelocityCur(:,iEdge) & + + weight1st * dt * normalVelocityTend(:,iEdge) + end do + ! Fine layers close to interface layers + do ie = 1, nEdgesInLTSRegion(1,3) + iEdge = edgesInLTSRegion(1,3,ie) + normalVelocityFirstStage(:,iEdge) = normalVelocityCur(:,iEdge) & + + weight1st * dt * normalVelocityTend(:,iEdge) + end do + + call mpas_timer_start("FB_LTS advance velocity solution") + call mpas_timer_start("FB_LTS advance solution") + ! END: Normal velocity stage 1 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! Halo update + call mpas_timer_start("FB_LTS prognostic halo update") + + call mpas_dmpar_field_halo_exch(domain, 'normalVelocity', timeLevel=3) + call mpas_dmpar_field_halo_exch(domain, 'layerThickness', timeLevel=3) + + call mpas_timer_stop("FB_LTS prognostic halo update") + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! BEGIN: Thickness stage 2 From 17b10d20ac33d062acb37b55c41386840cfe76fd Mon Sep 17 00:00:00 2001 From: Jeremy Lilly Date: Tue, 24 Oct 2023 09:04:14 -0700 Subject: [PATCH 11/34] Split prognostic halo update to do thickness and velocity seperately --- .../mpas_ocn_time_integration_fblts.F | 20 ++++++++++++++----- 1 file changed, 15 insertions(+), 5 deletions(-) diff --git a/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_fblts.F b/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_fblts.F index 5d6ccb484dc9..5bc421196cec 100644 --- a/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_fblts.F +++ b/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_fblts.F @@ -448,10 +448,19 @@ subroutine ocn_time_integrator_fblts(domain, dt)!{{{ call mpas_timer_stop("FB_LTS advance thickness solution") call mpas_timer_stop("FB_LTS advance solution") + ! Halo update + call mpas_timer_start("FB_LTS prognostic halo update") + call mpas_timer_start("FB_LTS thickness halo update") + + call mpas_dmpar_field_halo_exch(domain, 'layerThickness', timeLevel=3) + + call mpas_timer_stop("FB_LTS thickness halo update") + call mpas_timer_stop("FB_LTS prognostic halo update") + ! END: Thickness stage 1 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! BEGIN: Normal velocity stage 1 ! Compute the first stage of the normal velocity on fine*, @@ -498,17 +507,18 @@ subroutine ocn_time_integrator_fblts(domain, dt)!{{{ call mpas_timer_start("FB_LTS advance velocity solution") call mpas_timer_start("FB_LTS advance solution") - ! END: Normal velocity stage 1 - !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - ! Halo update call mpas_timer_start("FB_LTS prognostic halo update") + call mpas_timer_start("FB_LTS velocity halo update") call mpas_dmpar_field_halo_exch(domain, 'normalVelocity', timeLevel=3) - call mpas_dmpar_field_halo_exch(domain, 'layerThickness', timeLevel=3) + call mpas_timer_stop("FB_LTS velocity halo update") call mpas_timer_stop("FB_LTS prognostic halo update") + ! END: Normal velocity stage 1 + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! BEGIN: Thickness stage 2 From c5d6d8c67d73350510c9476cc5913019fbff1628 Mon Sep 17 00:00:00 2001 From: Jeremy Lilly Date: Fri, 27 Oct 2023 10:55:06 -0700 Subject: [PATCH 12/34] Implement FB averaging in vel_tend and write coarse advancement for stage 2, 3 --- .../mpas_ocn_time_integration_fblts.F | 299 ++++++++++++++++-- 1 file changed, 270 insertions(+), 29 deletions(-) diff --git a/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_fblts.F b/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_fblts.F index 5bc421196cec..49ad9c39c594 100644 --- a/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_fblts.F +++ b/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_fblts.F @@ -32,11 +32,11 @@ module ocn_time_integration_fblts use ocn_constants use ocn_tendency - use ocn_diagnostics + !use ocn_diagnostics use ocn_mesh use ocn_vmix use ocn_config - use ocn_diagnostics_variables + !use ocn_diagnostics_variables use ocn_equation_of_state use ocn_time_average_coupled use ocn_time_varying_forcing @@ -389,16 +389,6 @@ subroutine ocn_time_integrator_fblts(domain, dt)!{{{ ! BEGIN COARSE ADVANCEMENT !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - ! Update halos for diagnostic variables - call mpas_timer_start("FB_LTS diagnostic halo update") - - call mpas_dmpar_field_halo_exch(domain, 'normalizedRelativeVorticityEdge') - if (config_mom_del4 > 0.0_RKIND) then - call mpas_dmpar_field_halo_exch(domain, 'divergence') - call mpas_dmpar_field_halo_exch(domain, 'relativeVorticity') - end if - - call mpas_timer_stop("FB_LTS diagnostic halo update") !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! BEGIN: Thickness stage 1 @@ -419,7 +409,7 @@ subroutine ocn_time_integrator_fblts(domain, dt)!{{{ call mpas_timer_start("FB_LTS advance thickness solution") ! Interface layers - do iRegion =1,nRegions + do iRegion = 1,nRegions do ic = 1, nCellsInLTSRegion(iRegion,2) iCell = cellsInLTSRegion(iRegion,2,ic) do k = 1, maxLevelCell(iCell) @@ -466,11 +456,10 @@ subroutine ocn_time_integrator_fblts(domain, dt)!{{{ ! Compute the first stage of the normal velocity on fine*, ! interface 1, interface 2, and coarse - ! Compute fast tendencies for thickness + ! Compute fast tendencies for normal velocity call mpas_timer_start("FB_LTS compute fast tendencies") call mpas_timer_start("FB_LTS compute fast velocity tendencies") - ! NEED TO CALCULATE h^* call ocn_lts_vel_tend(statePool, LTSPool, tendPool, 1, 0, 1, 1, 1) call mpas_timer_stop("FB_LTS compute fast velocity tendencies") @@ -484,7 +473,7 @@ subroutine ocn_time_integrator_fblts(domain, dt)!{{{ + normalVelocityTendSlow(:,:) ! Interface regions - do iRegion =1,nRegions + do iRegion = 1,nRegions do ie = 1, nEdgesInLTSRegion(iRegion,2) iEdge = edgesInLTSRegion(iRegion,2,ie) normalVelocityFirstStage(:,iEdge) = normalVelocityCur(:,iEdge) & @@ -524,7 +513,59 @@ subroutine ocn_time_integrator_fblts(domain, dt)!{{{ ! BEGIN: Thickness stage 2 ! Compute the second stage of the thickness on fine*, interface 1, ! interface 2, and coarse + + ! Compute fast tendencies for thickness + call mpas_timer_start("FB_LTS compute fast tendencies") + call mpas_timer_start("FB_LTS compute fast thickness tendencies") + + call ocn_lts_thick_tend(statePool, LTSPool, tendPool, 3, 0, 1, 1, 1) + + call mpas_timer_stop("FB_LTS compute fast thickness tendencies") + call mpas_timer_stop("FB_LTS compute fast tendencies") + + ! Advance thickness solution with first stage + call mpas_timer_start("FB_LTS advance solution") + call mpas_timer_start("FB_LTS advance thickness solution") + ! Interface layers + do iRegion = 1,nRegions + do ic = 1, nCellsInLTSRegion(iRegion,2) + iCell = cellsInLTSRegion(iRegion,2,ic) + do k = 1, maxLevelCell(iCell) + layerThicknessSecondStage(k,iCell) = layerThicknessCur(k,iCell) & + + weight2nd * dt * layerThicknessTend(k,iCell) + end do + end do + end do + ! Coarse + do ic = 1, nCellsInLTSRegion(2,1) + iCell = cellsInLTSRegion(2,1,ic) + do k = 1, maxLevelCell(iCell) + layerThicknessSecondStage(k,iCell) = layerThicknessCur(k,iCell) & + + weight2nd * dt * layerThicknessTend(k,iCell) + end do + end do + ! Fine layers close to interface layers + do ic = 1, nCellsInLTSRegion(1,3) + iCell = cellsInLTSRegion(1,3,ic) + do k = 1, maxLevelCell(iCell) + layerThicknessSecondStage(k,iCell) = layerThicknessCur(k,iCell) & + + weight2nd * dt * layerThicknessTend(k,iCell) + end do + end do + + call mpas_timer_stop("FB_LTS advance thickness solution") + call mpas_timer_stop("FB_LTS advance solution") + + ! Halo update + call mpas_timer_start("FB_LTS prognostic halo update") + call mpas_timer_start("FB_LTS thickness halo update") + + call mpas_dmpar_field_halo_exch(domain, 'layerThickness', timeLevel=4) + + call mpas_timer_stop("FB_LTS thickness halo update") + call mpas_timer_stop("FB_LTS prognostic halo update") + ! END: Thickness stage 2 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! @@ -534,6 +575,55 @@ subroutine ocn_time_integrator_fblts(domain, dt)!{{{ ! Compute the second stage of the normal velocity on fine*, ! interface 1, interface 2, and coarse + ! Compute fast tendencies for normal velocity + call mpas_timer_start("FB_LTS compute fast tendencies") + call mpas_timer_start("FB_LTS compute fast velocity tendencies") + + call ocn_lts_vel_tend(statePool, LTSPool, tendPool, 2, 0, 1, 1, 1) + + call mpas_timer_stop("FB_LTS compute fast velocity tendencies") + call mpas_timer_stop("FB_LTS compute fast tendencies") + + ! Advance normal velocity solution with first stage + call mpas_timer_start("FB_LTS advance solution") + call mpas_timer_start("FB_LTS advance velocity solution") + + normalVelocityTend(:,:) = normalVelocityTend(:,:) & + + normalVelocityTendSlow(:,:) + + ! Interface regions + do iRegion = 1,nRegions + do ie = 1, nEdgesInLTSRegion(iRegion,2) + iEdge = edgesInLTSRegion(iRegion,2,ie) + normalVelocitySecondStage(:,iEdge) = normalVelocityCur(:,iEdge) & + + weight2nd * dt * normalVelocityTend(:,iEdge) + end do + end do + ! Coarse + do ie = 1, nEdgesInLTSRegion(2,1) + iEdge = edgesInLTSRegion(2,1,ie) + normalVelocitySecondStage(:,iEdge) = normalVelocityCur(:,iEdge) & + + weight2nd * dt * normalVelocityTend(:,iEdge) + end do + ! Fine layers close to interface layers + do ie = 1, nEdgesInLTSRegion(1,3) + iEdge = edgesInLTSRegion(1,3,ie) + normalVelocitySecondStage(:,iEdge) = normalVelocityCur(:,iEdge) & + + weight2nd * dt * normalVelocityTend(:,iEdge) + end do + + call mpas_timer_start("FB_LTS advance velocity solution") + call mpas_timer_start("FB_LTS advance solution") + + ! Halo update + call mpas_timer_start("FB_LTS prognostic halo update") + call mpas_timer_start("FB_LTS velocity halo update") + + call mpas_dmpar_field_halo_exch(domain, 'normalVelocity', timeLevel=4) + + call mpas_timer_stop("FB_LTS velocity halo update") + call mpas_timer_stop("FB_LTS prognostic halo update") + ! END: Normal velocity stage 2 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! @@ -542,7 +632,59 @@ subroutine ocn_time_integrator_fblts(domain, dt)!{{{ ! BEGIN: Thickness stage 3 ! Compute the third stage of the thickness on fine*, interface 1, ! interface 2, and coarse + + ! Compute fast tendencies for thickness + call mpas_timer_start("FB_LTS compute fast tendencies") + call mpas_timer_start("FB_LTS compute fast thickness tendencies") + + call ocn_lts_thick_tend(statePool, LTSPool, tendPool, 4, 0, 1, 1, 1) + + call mpas_timer_stop("FB_LTS compute fast thickness tendencies") + call mpas_timer_stop("FB_LTS compute fast tendencies") + + ! Advance thickness solution with first stage + call mpas_timer_start("FB_LTS advance solution") + call mpas_timer_start("FB_LTS advance thickness solution") + + ! Interface layers + do iRegion = 1,nRegions + do ic = 1, nCellsInLTSRegion(iRegion,2) + iCell = cellsInLTSRegion(iRegion,2,ic) + do k = 1, maxLevelCell(iCell) + layerThicknessNew(k,iCell) = layerThicknessCur(k,iCell) & + + dt * layerThicknessTend(k,iCell) + end do + end do + end do + ! Coarse + do ic = 1, nCellsInLTSRegion(2,1) + iCell = cellsInLTSRegion(2,1,ic) + do k = 1, maxLevelCell(iCell) + layerThicknessNew(k,iCell) = layerThicknessCur(k,iCell) & + + dt * layerThicknessTend(k,iCell) + end do + end do + ! Fine layers close to interface layers + do ic = 1, nCellsInLTSRegion(1,3) + iCell = cellsInLTSRegion(1,3,ic) + do k = 1, maxLevelCell(iCell) + layerThicknessNew(k,iCell) = layerThicknessCur(k,iCell) & + + dt * layerThicknessTend(k,iCell) + end do + end do + + call mpas_timer_stop("FB_LTS advance thickness solution") + call mpas_timer_stop("FB_LTS advance solution") + + ! Halo update + call mpas_timer_start("FB_LTS prognostic halo update") + call mpas_timer_start("FB_LTS thickness halo update") + call mpas_dmpar_field_halo_exch(domain, 'layerThickness', timeLevel=2) + + call mpas_timer_stop("FB_LTS thickness halo update") + call mpas_timer_stop("FB_LTS prognostic halo update") + ! END: Thickness stage 3 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! @@ -551,7 +693,58 @@ subroutine ocn_time_integrator_fblts(domain, dt)!{{{ ! BEGIN: Normal velocity stage 3 ! Compute the third stage of the normal velocity on, interface 1, ! and coarse + ! For now we compute on interface 2 as well, but this can/should + ! be removed later + + ! Compute fast tendencies for normal velocity + call mpas_timer_start("FB_LTS compute fast tendencies") + call mpas_timer_start("FB_LTS compute fast velocity tendencies") + + call ocn_lts_vel_tend(statePool, LTSPool, tendPool, 3, 0, 0, 1, 1) + + call mpas_timer_stop("FB_LTS compute fast velocity tendencies") + call mpas_timer_stop("FB_LTS compute fast tendencies") + + ! Advance normal velocity solution with first stage + call mpas_timer_start("FB_LTS advance solution") + call mpas_timer_start("FB_LTS advance velocity solution") + + normalVelocityTend(:,:) = normalVelocityTend(:,:) & + + normalVelocityTendSlow(:,:) + + ! Interface regions + do iRegion = 1,nRegions + do ie = 1, nEdgesInLTSRegion(iRegion,2) + iEdge = edgesInLTSRegion(iRegion,2,ie) + normalVelocityNew(:,iEdge) = normalVelocityCur(:,iEdge) & + + dt * normalVelocityTend(:,iEdge) + end do + end do + ! Coarse + do ie = 1, nEdgesInLTSRegion(2,1) + iEdge = edgesInLTSRegion(2,1,ie) + normalVelocityNew(:,iEdge) = normalVelocityCur(:,iEdge) & + + dt * normalVelocityTend(:,iEdge) + end do + ! Fine layers close to interface layers + !do ie = 1, nEdgesInLTSRegion(1,3) + ! iEdge = edgesInLTSRegion(1,3,ie) + ! normalVelocityNew(:,iEdge) = normalVelocityCur(:,iEdge) & + ! + dt * normalVelocityTend(:,iEdge) + !end do + + call mpas_timer_start("FB_LTS advance velocity solution") + call mpas_timer_start("FB_LTS advance solution") + + ! Halo update + call mpas_timer_start("FB_LTS prognostic halo update") + call mpas_timer_start("FB_LTS velocity halo update") + + call mpas_dmpar_field_halo_exch(domain, 'normalVelocity', timeLevel=2) + call mpas_timer_stop("FB_LTS velocity halo update") + call mpas_timer_stop("FB_LTS prognostic halo update") + ! END: Normal velocity stage 3 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! @@ -896,6 +1089,7 @@ subroutine ocn_lts_thick_tend(statePool, LTSPool, tendPool, timeLevelIn, compute call mpas_pool_get_array(statePool, 'layerThickness', layerThickness, timeLevelIn) ! calculate layerThickEdgeFlux + ! add upwind formulation, ask Giacomo about branch do iEdge = 1, nEdgesAll kmin = minLevelEdgeBot(iEdge) kmax = maxLevelEdgeTop(iEdge) @@ -986,14 +1180,14 @@ end subroutine ocn_lts_thick_tend!}}} ! ocn_lts_vel_tend ! !> \brief Calculate velocity tendencies for FB_LTS -!> \author Giacomo Capodaglio +!> \author Jeremy Lilly !> \date October 2023 !> \details !> This routine calculates the tendencies on different LTS regions ! !----------------------------------------------------------------------- - subroutine ocn_lts_vel_tend(statePool, LTSPool, tendPool, timeLevelIn, computeOnFineBig, & + subroutine ocn_lts_vel_tend(statePool, LTSPool, tendPool, stage, computeOnFineBig, & computeOnFineSmall, computeOnCoarse, computeOnInterface)!{{{ !----------------------------------------------------------------- @@ -1001,7 +1195,7 @@ subroutine ocn_lts_vel_tend(statePool, LTSPool, tendPool, timeLevelIn, computeOn !----------------------------------------------------------------- integer, intent(in) :: & - timeLevelIn, & + stage, & computeOnFineBig, & computeOnFineSmall, & computeOnCoarse, & @@ -1036,11 +1230,16 @@ subroutine ocn_lts_vel_tend(statePool, LTSPool, tendPool, timeLevelIn, computeOn real (kind=RKIND), dimension(:,:), pointer :: & normalVelocityTend, & normalVelocity, & - layerThickness + layerThickness, & + layerThicknessCur, & + layerThicknessNew, & + layerThicknessFirstStage, & + layerThicknessSecondStage integer :: & iEdge, cell1, cell2, k, ie, iRegion, nRegions, & - ic, i, iCell, kmin, kmax + ic, i, iCell, kmin, kmax, timeLevelVelocity + real (kind=RKIND) :: & invdcEdge, betaSelfAttrLoad, flux, ssh_sal_on @@ -1051,14 +1250,57 @@ subroutine ocn_lts_vel_tend(statePool, LTSPool, tendPool, timeLevelIn, computeOn betaSelfAttrLoad = config_self_attraction_and_loading_beta nRegions = 2 + ! Depending on which RK stage, calculate the forward-backward + ! average of the layer thickness to use in the tendency + ! calculations + ! Stage 1: h = weight1*firstStage + (1-weight1)*current + ! Stage 2: h = weight2*secondStage + (1-weight2)*current + ! Stage 3: h = weight3*new + (1-2*weight3)*secondStage + weight3*current + ! Also, set the desired timeLevel for normalVelocity + if (stage == 1) then + timeLevelVelocity = 1 + + call mpas_pool_get_array(statePool, 'layerThickness', & + layerThicknessCur, 1) + call mpas_pool_get_array(statePool, 'layerThickness', & + layerThicknessFirstStage, 3) + + layerThickness(:,:) = config_fb_weight_1 * layerThicknessFirstStage(:,:) & + + (1 - config_fb_weight_1) * layerThicknessCur(:,:) + else if (stage == 2) then + timeLevelVelocity = 3 + + call mpas_pool_get_array(statePool, 'layerThickness', & + layerThicknessCur, 1) + call mpas_pool_get_array(statePool, 'layerThickness', & + layerThicknessSecondStage, 4) + + layerThickness(:,:) = config_fb_weight_2 * layerThicknessSecondStage(:,:) & + + (1 - config_fb_weight_2) * layerThicknessCur(:,:) + else if (stage == 3) then + timeLevelVelocity = 4 + + call mpas_pool_get_array(statePool, 'layerThickness', & + layerThicknessCur, 1) + call mpas_pool_get_array(statePool, 'layerThickness', & + layerThicknessSecondStage, 4) + call mpas_pool_get_array(statePool, 'layerThickness', & + layerThicknessNew, 2) + + layerThickness(:,:) = config_fb_weight_3 * layerThicknessNew(:,:) & + + (1 - 2*config_fb_weight_3) * layerThicknessSecondStage(:,:) & + + config_fb_weight_3 * layerThicknessCur(:,:) + end if + call mpas_pool_get_array(LTSPool, 'edgesInLTSRegion', edgesInLTSRegion) call mpas_pool_get_array(LTSPool, 'nEdgesInLTSRegion', nEdgesInLTSRegion) call mpas_pool_get_array(tendPool, 'normalVelocity', normalVelocityTend) - call mpas_pool_get_array(statePool, 'normalVelocity', normalVelocity, timeLevelIn) - call mpas_pool_get_array(statePool, 'layerThickness', layerThickness, timeLevelIn) - call mpas_pool_get_array(statePool, 'ssh', ssh, timeLevelIn) + call mpas_pool_get_array(statePool, 'normalVelocity', & + normalVelocity, timeLevelVelocity) + call mpas_pool_get_array(statePool, 'ssh', & + ssh, timeLevelVelocity) if (config_use_self_attraction_loading) then ssh_sal_on = 1.0_RKIND @@ -1079,10 +1321,9 @@ subroutine ocn_lts_vel_tend(statePool, LTSPool, tendPool, timeLevelIn, computeOn normalVelocityTend(:,:) = 0.0_RKIND - ! interface + ! Interface regions if (computeOnInterface == 1) then do iRegion = 1, nRegions - ! velocity tendency do ie = 1, nEdgesInLTSRegion(iRegion,2) iEdge = edgesInLTSRegion(iRegion,2,ie) cell1 = cellsOnEdge(1,iEdge) @@ -1101,8 +1342,8 @@ subroutine ocn_lts_vel_tend(statePool, LTSPool, tendPool, timeLevelIn, computeOn end do end if + ! Fine in the interior, away from interface 1 if (computeOnFineBig == 1) then - ! velocity tendency do ie = 1, nEdgesInLTSRegion(1,1) iEdge = edgesInLTSRegion(1,1,ie) cell1 = cellsOnEdge(1,iEdge) @@ -1120,8 +1361,8 @@ subroutine ocn_lts_vel_tend(statePool, LTSPool, tendPool, timeLevelIn, computeOn end do end if + ! Fine adjacent to interface 1 if (computeOnFineSmall == 1) then - ! velocity tendency do ie = 1, nEdgesInLTSRegion(1,3) iEdge = edgesInLTSRegion(1,3,ie) cell1 = cellsOnEdge(1,iEdge) @@ -1139,8 +1380,8 @@ subroutine ocn_lts_vel_tend(statePool, LTSPool, tendPool, timeLevelIn, computeOn end do end if + ! Coarse if (computeOnCoarse == 1) then - ! velocity tendency do ie = 1, nEdgesInLTSRegion(2,1) iEdge = edgesInLTSRegion(2,1,ie) cell1 = cellsOnEdge(1,iEdge) From d63c768e2f4ba8e1fbd4af7c9f69b2b2bb875de2 Mon Sep 17 00:00:00 2001 From: Jeremy Lilly Date: Mon, 30 Oct 2023 11:14:25 -0700 Subject: [PATCH 13/34] Rewrite lts_thick_tend in preparation for calculating tends during subcycling --- .../mpas_ocn_time_integration_fblts.F | 104 ++++++++++-------- 1 file changed, 60 insertions(+), 44 deletions(-) diff --git a/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_fblts.F b/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_fblts.F index 49ad9c39c594..1952b37aa942 100644 --- a/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_fblts.F +++ b/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_fblts.F @@ -399,7 +399,9 @@ subroutine ocn_time_integrator_fblts(domain, dt)!{{{ call mpas_timer_start("FB_LTS compute fast tendencies") call mpas_timer_start("FB_LTS compute fast thickness tendencies") - call ocn_lts_thick_tend(statePool, LTSPool, tendPool, 1, 0, 1, 1, 1) + call ocn_lts_thick_tend(LTSPool, tendPool, & + normalVelocityCur, layerThicknessCur, & + 0, 1, 1, 1, 1) call mpas_timer_stop("FB_LTS compute fast thickness tendencies") call mpas_timer_stop("FB_LTS compute fast tendencies") @@ -518,7 +520,9 @@ subroutine ocn_time_integrator_fblts(domain, dt)!{{{ call mpas_timer_start("FB_LTS compute fast tendencies") call mpas_timer_start("FB_LTS compute fast thickness tendencies") - call ocn_lts_thick_tend(statePool, LTSPool, tendPool, 3, 0, 1, 1, 1) + call ocn_lts_thick_tend(LTSPool, tendPool, & + normalVelocityFirstStage, layerThicknessFirstStage, & + 0, 1, 1, 1, 1) call mpas_timer_stop("FB_LTS compute fast thickness tendencies") call mpas_timer_stop("FB_LTS compute fast tendencies") @@ -637,7 +641,9 @@ subroutine ocn_time_integrator_fblts(domain, dt)!{{{ call mpas_timer_start("FB_LTS compute fast tendencies") call mpas_timer_start("FB_LTS compute fast thickness tendencies") - call ocn_lts_thick_tend(statePool, LTSPool, tendPool, 4, 0, 1, 1, 1) + call ocn_lts_thick_tend(LTSPool, tendPool, & + normalVelocitySecondStage, layerThicknessSecondStage, & + 0, 1, 1, 1, 1) call mpas_timer_stop("FB_LTS compute fast thickness tendencies") call mpas_timer_stop("FB_LTS compute fast tendencies") @@ -1019,29 +1025,33 @@ end subroutine ocn_time_integration_fblts_init!}}} ! ocn_lts_thick_tend ! !> \brief Calculate thickness tendencies for FB_LTS -!> \author Giacomo Capodaglio +!> \author Jeremy Lilly !> \date October 2023 !> \details !> This routine calculates the tendencies on different LTS regions ! !----------------------------------------------------------------------- - subroutine ocn_lts_thick_tend(statePool, LTSPool, tendPool, timeLevelIn, computeOnFineBig, & - computeOnFineSmall, computeOnCoarse, computeOnInterface)!{{{ + subroutine ocn_lts_thick_tend(LTSPool, tendPool, & + normalVelocity, layerThickness, & + computeOnFineInterior, computeOnFineInterfaceAdjacent, & + computeOnInterfaceOne, computeOnInterfaceTwo, & + computeOnCoarse)!{{{ !----------------------------------------------------------------- ! Input variables !----------------------------------------------------------------- integer, intent(in) :: & - timeLevelIn, & - computeOnFineBig, & - computeOnFineSmall, & - computeOnCoarse, & - computeOnInterface + computeOnFineInterior, & + computeOnFineInterfaceAdjacent, & + computeOnInterfaceOne, & + computeOnInterfaceTwo, & + computeOnCoarse - type (mpas_pool_type), intent(in) :: & - statePool !< Input: state variables + real (kind=RKIND), dimension(:,:), pointer :: & + layerThickness, & + normalVelocity type (mpas_pool_type), intent(in) :: & LTSPool !< Input: LTS data @@ -1064,12 +1074,10 @@ subroutine ocn_lts_thick_tend(statePool, LTSPool, tendPool, timeLevelIn, compute cellsInLTSRegion real (kind=RKIND), dimension(:,:), pointer :: & - layerThicknessTend, & - normalVelocity, & - layerThickness + layerThicknessTend integer :: & - iEdge, cell1, cell2, k, ie, iRegion, nRegions, & + iEdge, cell1, cell2, k, ie, & ic, i, iCell, kmin, kmax real (kind=RKIND) :: & invdcEdge, flux @@ -1078,16 +1086,11 @@ subroutine ocn_lts_thick_tend(statePool, LTSPool, tendPool, timeLevelIn, compute ! Begin routine !----------------------------------------------------------------- - nRegions = 2 - call mpas_pool_get_array(LTSPool, 'cellsInLTSRegion', cellsInLTSRegion) call mpas_pool_get_array(LTSPool, 'nCellsInLTSRegion', nCellsInLTSRegion) call mpas_pool_get_array(tendPool, 'layerThickness', layerThicknessTend) - call mpas_pool_get_array(statePool, 'normalVelocity', normalVelocity, timeLevelIn) - call mpas_pool_get_array(statePool, 'layerThickness', layerThickness, timeLevelIn) - ! calculate layerThickEdgeFlux ! add upwind formulation, ask Giacomo about branch do iEdge = 1, nEdgesAll @@ -1110,25 +1113,8 @@ subroutine ocn_lts_thick_tend(statePool, LTSPool, tendPool, timeLevelIn, compute layerThicknessTend(:, :) = 0.0_RKIND - if (computeOnInterface == 1) then - do iRegion = 1, nRegions - ! thickness tendency - do ic = 1, nCellsInLTSRegion(iRegion,2) - iCell = cellsInLTSRegion(iRegion,2,ic) - do i = 1, nEdgesOnCell(iCell) - iEdge = edgesOnCell(i, iCell) - do k = minLevelEdgeBot(iEdge), maxLevelEdgeTop(iEdge) - flux = normalVelocity(k, iEdge) * dvEdge(iEdge) * layerThickEdgeFlux(k, iEdge) - layerThicknessTend(k, iCell) = layerThicknessTend(k, iCell) & - + edgeSignOnCell(i, iCell) * flux * invAreaCell(iCell) - end do - end do - end do - end do - end if - - if (computeOnFineBig == 1) then - ! thickness tendency + ! Fine in the interior, away from interface 1 + if (computeOnFineInterior == 1) then do ic = 1, nCellsInLTSRegion(1,1) iCell = cellsInLTSRegion(1,1,ic) do i = 1, nEdgesOnCell(iCell) @@ -1142,8 +1128,8 @@ subroutine ocn_lts_thick_tend(statePool, LTSPool, tendPool, timeLevelIn, compute end do end if - if (computeOnFineSmall == 1) then - ! thickness tendency + ! Fine adjacent to interface 1 + if (computeOnFineInterfaceAdjacent == 1) then do ic = 1, nCellsInLTSRegion(1,3) iCell = cellsInLTSRegion(1,3,ic) do i = 1, nEdgesOnCell(iCell) @@ -1157,8 +1143,38 @@ subroutine ocn_lts_thick_tend(statePool, LTSPool, tendPool, timeLevelIn, compute end do end if + ! Interface 1 + if (computeOnInterfaceOne == 1) then + do ic = 1, nCellsInLTSRegion(1,2) + iCell = cellsInLTSRegion(1,2,ic) + do i = 1, nEdgesOnCell(iCell) + iEdge = edgesOnCell(i, iCell) + do k = minLevelEdgeBot(iEdge), maxLevelEdgeTop(iEdge) + flux = normalVelocity(k, iEdge) * dvEdge(iEdge) * layerThickEdgeFlux(k, iEdge) + layerThicknessTend(k, iCell) = layerThicknessTend(k, iCell) & + + edgeSignOnCell(i, iCell) * flux * invAreaCell(iCell) + end do + end do + end do + end if + + ! Interface 2 + if (computeOnInterfaceTwo == 1) then + do ic = 1, nCellsInLTSRegion(2,2) + iCell = cellsInLTSRegion(2,2,ic) + do i = 1, nEdgesOnCell(iCell) + iEdge = edgesOnCell(i, iCell) + do k = minLevelEdgeBot(iEdge), maxLevelEdgeTop(iEdge) + flux = normalVelocity(k, iEdge) * dvEdge(iEdge) * layerThickEdgeFlux(k, iEdge) + layerThicknessTend(k, iCell) = layerThicknessTend(k, iCell) & + + edgeSignOnCell(i, iCell) * flux * invAreaCell(iCell) + end do + end do + end do + end if + + ! Coarse if (computeOnCoarse == 1) then - ! thickness tendency do ic = 1, nCellsInLTSRegion(2,1) iCell = cellsInLTSRegion(2,1,ic) do i = 1, nEdgesOnCell(iCell) From 5d2b872666194c54b46a3bb234a18e2427d6fcc6 Mon Sep 17 00:00:00 2001 From: Jeremy Lilly Date: Mon, 30 Oct 2023 15:50:02 -0700 Subject: [PATCH 14/34] Rewite lts_vel_tend similary to lts_thick_tend and put FB average calculations outside the tend routine --- .../mpas_ocn_time_integration_fblts.F | 188 +++++++++--------- 1 file changed, 93 insertions(+), 95 deletions(-) diff --git a/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_fblts.F b/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_fblts.F index 1952b37aa942..dba418454cb2 100644 --- a/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_fblts.F +++ b/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_fblts.F @@ -148,6 +148,11 @@ subroutine ocn_time_integrator_fblts(domain, dt)!{{{ layerThicknessNew, & ! layer thickness at time n+1 layerThicknessFirstStage, & ! layer thickness at first stage of LTS layerThicknessSecondStage ! layer thickness at second stage of LTS + + ! Extra pointers to store FB averages and interpolations of data stored in the state pool + real (kind=RKIND), dimension(:,:), pointer :: & + normalVelocityForTend, & + layerThicknessForTend ! LTS objects real (kind=RKIND) :: & @@ -458,11 +463,19 @@ subroutine ocn_time_integrator_fblts(domain, dt)!{{{ ! Compute the first stage of the normal velocity on fine*, ! interface 1, interface 2, and coarse + ! Perform forward-backward average of thickness for stage one of + ! FB-RK(3,2): + ! h = weight1*stage1 + (1-weight1)*current + layerThicknessForTend(:,:) = config_fb_weight_1 * layerThicknessFirstStage(:,:) & + + (1 - config_fb_weight_1) * layerThicknessCur(:,:) + ! Compute fast tendencies for normal velocity call mpas_timer_start("FB_LTS compute fast tendencies") call mpas_timer_start("FB_LTS compute fast velocity tendencies") - call ocn_lts_vel_tend(statePool, LTSPool, tendPool, 1, 0, 1, 1, 1) + call ocn_lts_vel_tend(LTSPool, tendPool, & + normalVelocityCur, layerThicknessForTend, & + 0, 1, 1, 1, 1) call mpas_timer_stop("FB_LTS compute fast velocity tendencies") call mpas_timer_stop("FB_LTS compute fast tendencies") @@ -579,11 +592,19 @@ subroutine ocn_time_integrator_fblts(domain, dt)!{{{ ! Compute the second stage of the normal velocity on fine*, ! interface 1, interface 2, and coarse + ! Perform forward-backward average of thickness for stage two of + ! FB-RK(3,2): + ! h = weight2*stage2 + (1-weight2)*current + layerThicknessForTend(:,:) = config_fb_weight_2 * layerThicknessSecondStage(:,:) & + + (1 - config_fb_weight_2) * layerThicknessCur(:,:) + ! Compute fast tendencies for normal velocity call mpas_timer_start("FB_LTS compute fast tendencies") call mpas_timer_start("FB_LTS compute fast velocity tendencies") - call ocn_lts_vel_tend(statePool, LTSPool, tendPool, 2, 0, 1, 1, 1) + call ocn_lts_vel_tend(LTSPool, tendPool, & + normalVelocityFirstStage, layerThicknessForTend, & + 0, 1, 1, 1, 1) call mpas_timer_stop("FB_LTS compute fast velocity tendencies") call mpas_timer_stop("FB_LTS compute fast tendencies") @@ -699,14 +720,21 @@ subroutine ocn_time_integrator_fblts(domain, dt)!{{{ ! BEGIN: Normal velocity stage 3 ! Compute the third stage of the normal velocity on, interface 1, ! and coarse - ! For now we compute on interface 2 as well, but this can/should - ! be removed later + ! Perform forward-backward average of thickness for stage one of + ! FB-RK(3,2): + ! h = weight3*new + (1-2*weight3)*secondStage + weight3*current + layerThicknessForTend(:,:) = config_fb_weight_3 * layerThicknessNew(:,:) & + + (1 - 2*config_fb_weight_3) * layerThicknessSecondStage(:,:) & + + config_fb_weight_3 * layerThicknessCur(:,:) + ! Compute fast tendencies for normal velocity call mpas_timer_start("FB_LTS compute fast tendencies") call mpas_timer_start("FB_LTS compute fast velocity tendencies") - call ocn_lts_vel_tend(statePool, LTSPool, tendPool, 3, 0, 0, 1, 1) + call ocn_lts_vel_tend(LTSPool, tendPool, & + normalVelocitySecondStage, layerThicknessForTend, & + 0, 0, 1, 0, 1) call mpas_timer_stop("FB_LTS compute fast velocity tendencies") call mpas_timer_stop("FB_LTS compute fast tendencies") @@ -1203,22 +1231,26 @@ end subroutine ocn_lts_thick_tend!}}} ! !----------------------------------------------------------------------- - subroutine ocn_lts_vel_tend(statePool, LTSPool, tendPool, stage, computeOnFineBig, & - computeOnFineSmall, computeOnCoarse, computeOnInterface)!{{{ + subroutine ocn_lts_vel_tend(LTSPool, tendPool, & + normalVelocity, layerThickness, & + computeOnFineInterior, computeOnFineInterfaceAdjacent, & + computeOnInterfaceOne, computeOnInterfaceTwo, & + computeOnCoarse)!{{{ !----------------------------------------------------------------- ! Input variables !----------------------------------------------------------------- integer, intent(in) :: & - stage, & - computeOnFineBig, & - computeOnFineSmall, & - computeOnCoarse, & - computeOnInterface + computeOnFineInterior, & + computeOnFineInterfaceAdjacent, & + computeOnInterfaceOne, & + computeOnInterfaceTwo, & + computeOnCoarse - type (mpas_pool_type), intent(in) :: & - statePool !< Input: state variables + real (kind=RKIND), dimension(:,:), pointer :: & + layerThickness, & + normalVelocity type (mpas_pool_type), intent(in) :: & LTSPool !< Input: LTS data @@ -1240,21 +1272,15 @@ subroutine ocn_lts_vel_tend(statePool, LTSPool, tendPool, stage, computeOnFineBi integer, dimension(:,:,:), pointer :: & edgesInLTSRegion - real (kind=RKIND), dimension(:), pointer :: & + real (kind=RKIND), dimension(nCellsAll) :: & ssh real (kind=RKIND), dimension(:,:), pointer :: & - normalVelocityTend, & - normalVelocity, & - layerThickness, & - layerThicknessCur, & - layerThicknessNew, & - layerThicknessFirstStage, & - layerThicknessSecondStage + normalVelocityTend integer :: & - iEdge, cell1, cell2, k, ie, iRegion, nRegions, & - ic, i, iCell, kmin, kmax, timeLevelVelocity + iEdge, cell1, cell2, k, ie, & + ic, i, iCell, kmin, kmax real (kind=RKIND) :: & invdcEdge, betaSelfAttrLoad, flux, ssh_sal_on @@ -1264,59 +1290,14 @@ subroutine ocn_lts_vel_tend(statePool, LTSPool, tendPool, stage, computeOnFineBi !----------------------------------------------------------------- betaSelfAttrLoad = config_self_attraction_and_loading_beta - nRegions = 2 - - ! Depending on which RK stage, calculate the forward-backward - ! average of the layer thickness to use in the tendency - ! calculations - ! Stage 1: h = weight1*firstStage + (1-weight1)*current - ! Stage 2: h = weight2*secondStage + (1-weight2)*current - ! Stage 3: h = weight3*new + (1-2*weight3)*secondStage + weight3*current - ! Also, set the desired timeLevel for normalVelocity - if (stage == 1) then - timeLevelVelocity = 1 - - call mpas_pool_get_array(statePool, 'layerThickness', & - layerThicknessCur, 1) - call mpas_pool_get_array(statePool, 'layerThickness', & - layerThicknessFirstStage, 3) - - layerThickness(:,:) = config_fb_weight_1 * layerThicknessFirstStage(:,:) & - + (1 - config_fb_weight_1) * layerThicknessCur(:,:) - else if (stage == 2) then - timeLevelVelocity = 3 - - call mpas_pool_get_array(statePool, 'layerThickness', & - layerThicknessCur, 1) - call mpas_pool_get_array(statePool, 'layerThickness', & - layerThicknessSecondStage, 4) - - layerThickness(:,:) = config_fb_weight_2 * layerThicknessSecondStage(:,:) & - + (1 - config_fb_weight_2) * layerThicknessCur(:,:) - else if (stage == 3) then - timeLevelVelocity = 4 - - call mpas_pool_get_array(statePool, 'layerThickness', & - layerThicknessCur, 1) - call mpas_pool_get_array(statePool, 'layerThickness', & - layerThicknessSecondStage, 4) - call mpas_pool_get_array(statePool, 'layerThickness', & - layerThicknessNew, 2) - - layerThickness(:,:) = config_fb_weight_3 * layerThicknessNew(:,:) & - + (1 - 2*config_fb_weight_3) * layerThicknessSecondStage(:,:) & - + config_fb_weight_3 * layerThicknessCur(:,:) - end if call mpas_pool_get_array(LTSPool, 'edgesInLTSRegion', edgesInLTSRegion) call mpas_pool_get_array(LTSPool, 'nEdgesInLTSRegion', nEdgesInLTSRegion) call mpas_pool_get_array(tendPool, 'normalVelocity', normalVelocityTend) - call mpas_pool_get_array(statePool, 'normalVelocity', & - normalVelocity, timeLevelVelocity) - call mpas_pool_get_array(statePool, 'ssh', & - ssh, timeLevelVelocity) + !call mpas_pool_get_array(statePool, 'ssh', & + ! ssh, timeLevelVelocity) if (config_use_self_attraction_loading) then ssh_sal_on = 1.0_RKIND @@ -1329,7 +1310,7 @@ subroutine ocn_lts_vel_tend(statePool, LTSPool, tendPool, stage, computeOnFineBi k = maxLevelCell(iCell) zTop(k:nVertLevels,iCell) = -bottomDepth(iCell) + layerThickness(k,iCell) do k = maxLevelCell(iCell)-1, minLevelCell(iCell), -1 - zTop(k,iCell) = zTop(k+1,iCell) + layerThickness(k ,iCell) + zTop(k,iCell) = zTop(k+1,iCell) + layerThickness(k,iCell) end do ! copy zTop(1,iCell) into sea-surface height array ssh(iCell) = zTop(minLevelCell(iCell),iCell) @@ -1337,29 +1318,8 @@ subroutine ocn_lts_vel_tend(statePool, LTSPool, tendPool, stage, computeOnFineBi normalVelocityTend(:,:) = 0.0_RKIND - ! Interface regions - if (computeOnInterface == 1) then - do iRegion = 1, nRegions - do ie = 1, nEdgesInLTSRegion(iRegion,2) - iEdge = edgesInLTSRegion(iRegion,2,ie) - cell1 = cellsOnEdge(1,iEdge) - cell2 = cellsOnEdge(2,iEdge) - invdcEdge = 1.0_RKIND / dcEdge(iEdge) - kMin = minLevelEdgeBot(iEdge) - kMax = maxLevelEdgeTop(iEdge) - do k=kMin,kMax - normalVelocityTend(k,iEdge) = normalVelocityTend(k,iEdge) & - - edgeMask(k,iEdge) * invdcEdge & - * ( gravity * ( (ssh(cell2) - ssh(cell1)) & - - (1.0_RKIND - ssh_sal_on) * betaSelfAttrLoad & - * (ssh(cell2) - ssh(cell1)) ) ) - end do - end do - end do - end if - ! Fine in the interior, away from interface 1 - if (computeOnFineBig == 1) then + if (computeOnFineInterior == 1) then do ie = 1, nEdgesInLTSRegion(1,1) iEdge = edgesInLTSRegion(1,1,ie) cell1 = cellsOnEdge(1,iEdge) @@ -1378,7 +1338,7 @@ subroutine ocn_lts_vel_tend(statePool, LTSPool, tendPool, stage, computeOnFineBi end if ! Fine adjacent to interface 1 - if (computeOnFineSmall == 1) then + if (computeOnFineInterfaceAdjacent == 1) then do ie = 1, nEdgesInLTSRegion(1,3) iEdge = edgesInLTSRegion(1,3,ie) cell1 = cellsOnEdge(1,iEdge) @@ -1396,6 +1356,44 @@ subroutine ocn_lts_vel_tend(statePool, LTSPool, tendPool, stage, computeOnFineBi end do end if + ! Interface 1 + if (computeOnInterfaceOne == 1) then + do ie = 1, nEdgesInLTSRegion(1,2) + iEdge = edgesInLTSRegion(1,2,ie) + cell1 = cellsOnEdge(1,iEdge) + cell2 = cellsOnEdge(2,iEdge) + invdcEdge = 1.0_RKIND / dcEdge(iEdge) + kMin = minLevelEdgeBot(iEdge) + kMax = maxLevelEdgeTop(iEdge) + do k=kMin,kMax + normalVelocityTend(k,iEdge) = normalVelocityTend(k,iEdge) & + - edgeMask(k,iEdge) * invdcEdge & + * ( gravity * ( (ssh(cell2) - ssh(cell1)) & + - (1.0_RKIND - ssh_sal_on) * betaSelfAttrLoad & + * (ssh(cell2) - ssh(cell1)) ) ) + end do + end do + end if + + ! Interface 2 + if (computeOnInterfaceTwo == 1) then + do ie = 1, nEdgesInLTSRegion(2,2) + iEdge = edgesInLTSRegion(2,2,ie) + cell1 = cellsOnEdge(1,iEdge) + cell2 = cellsOnEdge(2,iEdge) + invdcEdge = 1.0_RKIND / dcEdge(iEdge) + kMin = minLevelEdgeBot(iEdge) + kMax = maxLevelEdgeTop(iEdge) + do k=kMin,kMax + normalVelocityTend(k,iEdge) = normalVelocityTend(k,iEdge) & + - edgeMask(k,iEdge) * invdcEdge & + * ( gravity * ( (ssh(cell2) - ssh(cell1)) & + - (1.0_RKIND - ssh_sal_on) * betaSelfAttrLoad & + * (ssh(cell2) - ssh(cell1)) ) ) + end do + end do + end if + ! Coarse if (computeOnCoarse == 1) then do ie = 1, nEdgesInLTSRegion(2,1) From 7da1b00e7f6f9006ba014f13b75835b3588cdf25 Mon Sep 17 00:00:00 2001 From: Jeremy Lilly Date: Tue, 31 Oct 2023 14:36:51 -0700 Subject: [PATCH 15/34] Add a 5th time level to the statePool for storing averages of RK stage data, rearrage halo updates to only appear immediately before a tendency calculation --- .../mpas_ocn_time_integration_fblts.F | 188 ++++++++++++------ 1 file changed, 131 insertions(+), 57 deletions(-) diff --git a/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_fblts.F b/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_fblts.F index dba418454cb2..35d8e2cad7a6 100644 --- a/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_fblts.F +++ b/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_fblts.F @@ -144,16 +144,13 @@ subroutine ocn_time_integrator_fblts(domain, dt)!{{{ normalVelocityNew, & ! normal velocity at time n+1 normalVelocityFirstStage, & ! normal velocity at first stage of LTS normalVelocitySecondStage, & ! normal velocity at second stage of LTS + normalVelocityForTend, & ! extra variable to store averages of data for tends calculations layerThicknessCur, & ! layer thickness at time n layerThicknessNew, & ! layer thickness at time n+1 layerThicknessFirstStage, & ! layer thickness at first stage of LTS - layerThicknessSecondStage ! layer thickness at second stage of LTS + layerThicknessSecondStage, & ! layer thickness at second stage of LTS + layerThicknessForTend ! extra variable to store averages of data for tends calculations - ! Extra pointers to store FB averages and interpolations of data stored in the state pool - real (kind=RKIND), dimension(:,:), pointer :: & - normalVelocityForTend, & - layerThicknessForTend - ! LTS objects real (kind=RKIND) :: & dtFine, & ! fine dt, defined as dt / M @@ -216,6 +213,8 @@ subroutine ocn_time_integrator_fblts(domain, dt)!{{{ normalVelocityFirstStage, 3) call mpas_pool_get_array(statePool, 'normalVelocity', & normalVelocitySecondStage, 4) + call mpas_pool_get_array(statePool, 'normalVelocity', & + normalVelocityForTend, 5) call mpas_pool_get_array(statePool, 'layerThickness', & layerThicknessCur, 1) call mpas_pool_get_array(statePool, 'layerThickness', & @@ -224,6 +223,8 @@ subroutine ocn_time_integrator_fblts(domain, dt)!{{{ layerThicknessFirstStage, 3) call mpas_pool_get_array(statePool, 'layerThickness', & layerThicknessSecondStage, 4) + call mpas_pool_get_array(statePool, 'layerThickness', & + layerThicknessForTend, 5) ! Retrieve tendency variables call mpas_pool_get_array(tendPool, 'normalVelocity', normalVelocityTend) @@ -265,6 +266,10 @@ subroutine ocn_time_integrator_fblts(domain, dt)!{{{ call mpas_pool_get_array(tendSum3rdPool, 'layerThickness', & layerThicknessTendSum3rd) + ! not sure if we need this yet -- also need to dealloc if used! + !allocate(normalVelocityForTend(nVertLevels,nEdgesAll)) + !allocate(layerThicknessForTend(nVertLevels,nCellsAll)) + ! Init variables do iEdge = 1, nEdgesAll ! do we need this? do k = 1, maxLevelEdgeTop(iEdge) @@ -445,15 +450,6 @@ subroutine ocn_time_integrator_fblts(domain, dt)!{{{ call mpas_timer_stop("FB_LTS advance thickness solution") call mpas_timer_stop("FB_LTS advance solution") - ! Halo update - call mpas_timer_start("FB_LTS prognostic halo update") - call mpas_timer_start("FB_LTS thickness halo update") - - call mpas_dmpar_field_halo_exch(domain, 'layerThickness', timeLevel=3) - - call mpas_timer_stop("FB_LTS thickness halo update") - call mpas_timer_stop("FB_LTS prognostic halo update") - ! END: Thickness stage 1 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! @@ -468,6 +464,17 @@ subroutine ocn_time_integrator_fblts(domain, dt)!{{{ ! h = weight1*stage1 + (1-weight1)*current layerThicknessForTend(:,:) = config_fb_weight_1 * layerThicknessFirstStage(:,:) & + (1 - config_fb_weight_1) * layerThicknessCur(:,:) + + ! Halo update before tend calculation + ! Might be able to remove this? + call mpas_timer_start("FB_LTS prognostic halo update") + call mpas_timer_start("FB_LTS thickness halo update") + + ! Don't need to update normalVelocityCur, already updated + call mpas_dmpar_field_halo_exch(domain, 'layerThickness', timeLevel=5) + + call mpas_timer_stop("FB_LTS thickness halo update") + call mpas_timer_stop("FB_LTS prognostic halo update") ! Compute fast tendencies for normal velocity call mpas_timer_start("FB_LTS compute fast tendencies") @@ -511,15 +518,6 @@ subroutine ocn_time_integrator_fblts(domain, dt)!{{{ call mpas_timer_start("FB_LTS advance velocity solution") call mpas_timer_start("FB_LTS advance solution") - ! Halo update - call mpas_timer_start("FB_LTS prognostic halo update") - call mpas_timer_start("FB_LTS velocity halo update") - - call mpas_dmpar_field_halo_exch(domain, 'normalVelocity', timeLevel=3) - - call mpas_timer_stop("FB_LTS velocity halo update") - call mpas_timer_stop("FB_LTS prognostic halo update") - ! END: Normal velocity stage 1 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! @@ -529,6 +527,16 @@ subroutine ocn_time_integrator_fblts(domain, dt)!{{{ ! Compute the second stage of the thickness on fine*, interface 1, ! interface 2, and coarse + ! Halo update before tend calculation + call mpas_timer_start("FB_LTS prognostic halo update") + call mpas_timer_start("FB_LTS velocity halo update") + + call mpas_dmpar_field_halo_exch(domain, 'normalVelocity', timeLevel=3) + call mpas_dmpar_field_halo_exch(domain, 'layerThickness', timeLevel=3) + + call mpas_timer_stop("FB_LTS velocity halo update") + call mpas_timer_stop("FB_LTS prognostic halo update") + ! Compute fast tendencies for thickness call mpas_timer_start("FB_LTS compute fast tendencies") call mpas_timer_start("FB_LTS compute fast thickness tendencies") @@ -573,16 +581,7 @@ subroutine ocn_time_integrator_fblts(domain, dt)!{{{ call mpas_timer_stop("FB_LTS advance thickness solution") call mpas_timer_stop("FB_LTS advance solution") - - ! Halo update - call mpas_timer_start("FB_LTS prognostic halo update") - call mpas_timer_start("FB_LTS thickness halo update") - - call mpas_dmpar_field_halo_exch(domain, 'layerThickness', timeLevel=4) - call mpas_timer_stop("FB_LTS thickness halo update") - call mpas_timer_stop("FB_LTS prognostic halo update") - ! END: Thickness stage 2 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! @@ -598,6 +597,17 @@ subroutine ocn_time_integrator_fblts(domain, dt)!{{{ layerThicknessForTend(:,:) = config_fb_weight_2 * layerThicknessSecondStage(:,:) & + (1 - config_fb_weight_2) * layerThicknessCur(:,:) + ! Halo update before tend calculation + ! Might be able to remove this? + call mpas_timer_start("FB_LTS prognostic halo update") + call mpas_timer_start("FB_LTS thickness halo update") + + ! Don't need to update normalVelocityFirstStage, previously updated + call mpas_dmpar_field_halo_exch(domain, 'layerThickness', timeLevel=5) + + call mpas_timer_stop("FB_LTS thickness halo update") + call mpas_timer_stop("FB_LTS prognostic halo update") + ! Compute fast tendencies for normal velocity call mpas_timer_start("FB_LTS compute fast tendencies") call mpas_timer_start("FB_LTS compute fast velocity tendencies") @@ -640,15 +650,7 @@ subroutine ocn_time_integrator_fblts(domain, dt)!{{{ call mpas_timer_start("FB_LTS advance velocity solution") call mpas_timer_start("FB_LTS advance solution") - ! Halo update - call mpas_timer_start("FB_LTS prognostic halo update") - call mpas_timer_start("FB_LTS velocity halo update") - - call mpas_dmpar_field_halo_exch(domain, 'normalVelocity', timeLevel=4) - call mpas_timer_stop("FB_LTS velocity halo update") - call mpas_timer_stop("FB_LTS prognostic halo update") - ! END: Normal velocity stage 2 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! @@ -658,6 +660,16 @@ subroutine ocn_time_integrator_fblts(domain, dt)!{{{ ! Compute the third stage of the thickness on fine*, interface 1, ! interface 2, and coarse + ! Halo update before tends calculation + call mpas_timer_start("FB_LTS prognostic halo update") + call mpas_timer_start("FB_LTS velocity halo update") + + call mpas_dmpar_field_halo_exch(domain, 'normalVelocity', timeLevel=4) + call mpas_dmpar_field_halo_exch(domain, 'layerThickness', timeLevel=4) + + call mpas_timer_stop("FB_LTS velocity halo update") + call mpas_timer_stop("FB_LTS prognostic halo update") + ! Compute fast tendencies for thickness call mpas_timer_start("FB_LTS compute fast tendencies") call mpas_timer_start("FB_LTS compute fast thickness tendencies") @@ -703,15 +715,6 @@ subroutine ocn_time_integrator_fblts(domain, dt)!{{{ call mpas_timer_stop("FB_LTS advance thickness solution") call mpas_timer_stop("FB_LTS advance solution") - ! Halo update - call mpas_timer_start("FB_LTS prognostic halo update") - call mpas_timer_start("FB_LTS thickness halo update") - - call mpas_dmpar_field_halo_exch(domain, 'layerThickness', timeLevel=2) - - call mpas_timer_stop("FB_LTS thickness halo update") - call mpas_timer_stop("FB_LTS prognostic halo update") - ! END: Thickness stage 3 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! @@ -728,6 +731,17 @@ subroutine ocn_time_integrator_fblts(domain, dt)!{{{ + (1 - 2*config_fb_weight_3) * layerThicknessSecondStage(:,:) & + config_fb_weight_3 * layerThicknessCur(:,:) + ! Halo update before tend calculation + ! Might be able to remove this? + call mpas_timer_start("FB_LTS prognostic halo update") + call mpas_timer_start("FB_LTS thickness halo update") + + ! Don't need to update normalVelocitySecondStage, previously updated + call mpas_dmpar_field_halo_exch(domain, 'layerThickness', timeLevel=5) + + call mpas_timer_stop("FB_LTS thickness halo update") + call mpas_timer_stop("FB_LTS prognostic halo update") + ! Compute fast tendencies for normal velocity call mpas_timer_start("FB_LTS compute fast tendencies") call mpas_timer_start("FB_LTS compute fast velocity tendencies") @@ -760,17 +774,12 @@ subroutine ocn_time_integrator_fblts(domain, dt)!{{{ normalVelocityNew(:,iEdge) = normalVelocityCur(:,iEdge) & + dt * normalVelocityTend(:,iEdge) end do - ! Fine layers close to interface layers - !do ie = 1, nEdgesInLTSRegion(1,3) - ! iEdge = edgesInLTSRegion(1,3,ie) - ! normalVelocityNew(:,iEdge) = normalVelocityCur(:,iEdge) & - ! + dt * normalVelocityTend(:,iEdge) - !end do call mpas_timer_start("FB_LTS advance velocity solution") call mpas_timer_start("FB_LTS advance solution") ! Halo update + ! Might be able to remove this? call mpas_timer_start("FB_LTS prognostic halo update") call mpas_timer_start("FB_LTS velocity halo update") @@ -790,11 +799,76 @@ subroutine ocn_time_integrator_fblts(domain, dt)!{{{ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! BEGIN FINE ADVANCEMENT !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - do im = 1,M + do im = 0, M-1 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! BEGIN: Thickness stage 1 ! Compute the first stage of the thickness on fine + + ! Calculate predicted 'current' data at intermediate + ! time-level on interface 1: + ! (u,h) = (im/M)*new + (1- (im/M))*current + call mpas_timer_start('FB_LTS interface prediction') + + normalVelocityForTend(:,:) = normalVelocityCur(:,:) + do ie = 1, nEdgesInLTSRegion(1,2) + iEdge = edgesInLTSRegion(1,2,ie) + normalVelocityForTend(:,iEdge) = (im/M) * normalVelocityNew(:,iEdge) & + + (1 - (im/M)) * normalVelocityCur(:,iEdge) + end do + + layerThicknessForTend(:,:) = layerThicknessCur(:,:) + do ic = 1, nCellsInLTSRegion(1,2) + iCell = cellsInLTSRegion(1,2,ic) + layerThicknessForTend(:,iCell) = (im/M) * layerThicknessNew(:,iCell) & + + (1 - (im/M)) * layerThicknessCur(:,iCell) + end do + + call mpas_timer_stop('FB_LTS interface prediction') + + ! Compute fast tendencies for thickness + call mpas_timer_start("FB_LTS compute fast tendencies") + call mpas_timer_start("FB_LTS compute fast thickness tendencies") + + call ocn_lts_thick_tend(LTSPool, tendPool, & + normalVelocityForTend, layerThicknessForTend, & + 1, 1, 0, 0, 0) + + call mpas_timer_stop("FB_LTS compute fast thickness tendencies") + call mpas_timer_stop("FB_LTS compute fast tendencies") + + ! Advance thickness solution with first stage + call mpas_timer_start("FB_LTS advance solution") + call mpas_timer_start("FB_LTS advance thickness solution") + + ! Fine cells adjacent to interface 1 + do ic = 1, nCellsInLTSRegion(1,3) + iCell = cellsInLTSRegion(1,3,ic) + do k = 1, maxLevelCell(iCell) + layerThicknessFirstStage(k,iCell) = layerThicknessCur(k,iCell) & + + weight1st * dt * layerThicknessTend(k,iCell) + end do + end do + ! Fine in the interior, away from interface 1 + do ic = 1, nCellsInLTSRegion(1,1) + iCell = cellsInLTSRegion(1,1,ic) + do k = 1, maxLevelCell(iCell) + layerThicknessFirstStage(k,iCell) = layerThicknessCur(k,iCell) & + + weight1st * dt * layerThicknessTend(k,iCell) + end do + end do + + call mpas_timer_stop("FB_LTS advance thickness solution") + call mpas_timer_stop("FB_LTS advance solution") + + ! Halo update + !call mpas_timer_start("FB_LTS prognostic halo update") + !call mpas_timer_start("FB_LTS thickness halo update") + + !call mpas_dmpar_field_halo_exch(domain, 'layerThickness', timeLevel=1) + + !call mpas_timer_stop("FB_LTS thickness halo update") + !call mpas_timer_stop("FB_LTS prognostic halo update") ! END: Thickness stage 1 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! From 9f936e29f081b93b88cfa99c6f7f11ce62f6cc6d Mon Sep 17 00:00:00 2001 From: Jeremy Lilly Date: Tue, 31 Oct 2023 15:33:18 -0700 Subject: [PATCH 16/34] Implement the first stage for the fine subcycling --- .../mpas_ocn_time_integration_fblts.F | 107 ++++++++++++++---- 1 file changed, 87 insertions(+), 20 deletions(-) diff --git a/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_fblts.F b/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_fblts.F index 35d8e2cad7a6..6b0807e3d08e 100644 --- a/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_fblts.F +++ b/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_fblts.F @@ -515,8 +515,8 @@ subroutine ocn_time_integrator_fblts(domain, dt)!{{{ + weight1st * dt * normalVelocityTend(:,iEdge) end do - call mpas_timer_start("FB_LTS advance velocity solution") - call mpas_timer_start("FB_LTS advance solution") + call mpas_timer_stop("FB_LTS advance velocity solution") + call mpas_timer_stop("FB_LTS advance solution") ! END: Normal velocity stage 1 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! @@ -647,9 +647,8 @@ subroutine ocn_time_integrator_fblts(domain, dt)!{{{ + weight2nd * dt * normalVelocityTend(:,iEdge) end do - call mpas_timer_start("FB_LTS advance velocity solution") - call mpas_timer_start("FB_LTS advance solution") - + call mpas_timer_stop("FB_LTS advance velocity solution") + call mpas_timer_stop("FB_LTS advance solution") ! END: Normal velocity stage 2 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! @@ -724,7 +723,7 @@ subroutine ocn_time_integrator_fblts(domain, dt)!{{{ ! Compute the third stage of the normal velocity on, interface 1, ! and coarse - ! Perform forward-backward average of thickness for stage one of + ! Perform forward-backward average of thickness for stage three of ! FB-RK(3,2): ! h = weight3*new + (1-2*weight3)*secondStage + weight3*current layerThicknessForTend(:,:) = config_fb_weight_3 * layerThicknessNew(:,:) & @@ -775,8 +774,8 @@ subroutine ocn_time_integrator_fblts(domain, dt)!{{{ + dt * normalVelocityTend(:,iEdge) end do - call mpas_timer_start("FB_LTS advance velocity solution") - call mpas_timer_start("FB_LTS advance solution") + call mpas_timer_stop("FB_LTS advance velocity solution") + call mpas_timer_stop("FB_LTS advance solution") ! Halo update ! Might be able to remove this? @@ -826,6 +825,16 @@ subroutine ocn_time_integrator_fblts(domain, dt)!{{{ call mpas_timer_stop('FB_LTS interface prediction') + ! Halo update before tend calculation + call mpas_timer_start("FB_LTS prognostic halo update") + call mpas_timer_start("FB_LTS thickness halo update") + + call mpas_dmpar_field_halo_exch(domain, 'normalVelocity', timeLevel=5) + call mpas_dmpar_field_halo_exch(domain, 'layerThickness', timeLevel=5) + + call mpas_timer_stop("FB_LTS thickness halo update") + call mpas_timer_stop("FB_LTS prognostic halo update") + ! Compute fast tendencies for thickness call mpas_timer_start("FB_LTS compute fast tendencies") call mpas_timer_start("FB_LTS compute fast thickness tendencies") @@ -846,7 +855,7 @@ subroutine ocn_time_integrator_fblts(domain, dt)!{{{ iCell = cellsInLTSRegion(1,3,ic) do k = 1, maxLevelCell(iCell) layerThicknessFirstStage(k,iCell) = layerThicknessCur(k,iCell) & - + weight1st * dt * layerThicknessTend(k,iCell) + + weight1st * dtFine * layerThicknessTend(k,iCell) end do end do ! Fine in the interior, away from interface 1 @@ -854,22 +863,13 @@ subroutine ocn_time_integrator_fblts(domain, dt)!{{{ iCell = cellsInLTSRegion(1,1,ic) do k = 1, maxLevelCell(iCell) layerThicknessFirstStage(k,iCell) = layerThicknessCur(k,iCell) & - + weight1st * dt * layerThicknessTend(k,iCell) + + weight1st * dtFine * layerThicknessTend(k,iCell) end do end do call mpas_timer_stop("FB_LTS advance thickness solution") call mpas_timer_stop("FB_LTS advance solution") - ! Halo update - !call mpas_timer_start("FB_LTS prognostic halo update") - !call mpas_timer_start("FB_LTS thickness halo update") - - !call mpas_dmpar_field_halo_exch(domain, 'layerThickness', timeLevel=1) - - !call mpas_timer_stop("FB_LTS thickness halo update") - !call mpas_timer_stop("FB_LTS prognostic halo update") - ! END: Thickness stage 1 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! @@ -877,7 +877,74 @@ subroutine ocn_time_integrator_fblts(domain, dt)!{{{ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! BEGIN: Normal velocity stage 1 ! Compute the first stage of the normal velocity on fine - + + ! Perform forward-backward average of thickness for stage one of + ! FB-RK(3,2): + ! h = weight1*stage1 + (1-weight1)*current + layerThicknessForTend(:,:) = config_fb_weight_1 * layerThicknessFirstStage(:,:) & + + (1 - config_fb_weight_1) * layerThicknessCur(:,:) + + ! Calculate predicted FB averaged thickness data at + ! intermediate time-level on interface 1: + ! h = weight1*'firstStage' + (1-weight1)*'current' + ! = weight1*( (im/M)*new + (1/M)*firstStage + (1-(im+1)/M)*current ) + ! + (1-weight1)*( (im/M)*new + (1-im/M)*current ) + call mpas_timer_start('FB_LTS interface prediction') + + do ic = 1, nCellsInLTSRegion(1,2) + iCell = cellsInLTSRegion(1,2,ic) + layerThicknessForTend(:,iCell) = config_fb_weight_1 & + * ( (im/M) * layerThicknessNew(:,iCell) & + + (1/M) * layerThicknessFirstStage(:,iCell) & + + (1 - ((im+1)/M)) * layerThicknessCur(:,iCell) ) & + + (1 - config_fb_weight_1) & + * ( (im/M) * layerThicknessNew(:,iCell) & + + (1 - (im/M)) * layerThicknessCur(:,iCell) ) + end do + + call mpas_timer_stop('FB_LTS interface prediction') + + ! Halo update before tend calculation + call mpas_timer_start("FB_LTS prognostic halo update") + call mpas_timer_start("FB_LTS thickness halo update") + + ! Don't need to update normalVelocityForTend, previously updated + call mpas_dmpar_field_halo_exch(domain, 'layerThickness', timeLevel=5) + + call mpas_timer_stop("FB_LTS thickness halo update") + call mpas_timer_stop("FB_LTS prognostic halo update") + + ! Compute fast tendencies for normal velocity + call mpas_timer_start("FB_LTS compute fast tendencies") + call mpas_timer_start("FB_LTS compute fast velocity tendencies") + + call ocn_lts_vel_tend(LTSPool, tendPool, & + normalVelocityForTend, layerThicknessForTend, & + 1, 1, 0, 0, 0) + + call mpas_timer_stop("FB_LTS compute fast velocity tendencies") + call mpas_timer_stop("FB_LTS compute fast tendencies") + + ! Advance normal velocity solution with first stage + call mpas_timer_start("FB_LTS advance solution") + call mpas_timer_start("FB_LTS advance velocity solution") + + ! Fine cells adjacent to interface 1 + do ie = 1, nEdgesInLTSRegion(1,3) + iEdge = edgesInLTSRegion(1,3,ie) + normalVelocityFirstStage(:,iEdge) = normalVelocityCur(:,iEdge) & + + weight1st * dtFine * normalVelocityTend(:,iEdge) + end do + ! Fine in the interior, away from interface 1 + do ie = 1, nEdgesInLTSRegion(1,1) + iEdge = edgesInLTSRegion(1,1,ie) + normalVelocityFirstStage(:,iEdge) = normalVelocityCur(:,iEdge) & + + weight1st * dtFine * normalVelocityTend(:,iEdge) + end do + + call mpas_timer_stop("FB_LTS advance velocity solution") + call mpas_timer_stop("FB_LTS advance solution") + ! END: Normal velocity stage 1 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! From 1230a543ea15833bbc6705a7bf5aac1d45da0e3c Mon Sep 17 00:00:00 2001 From: Jeremy Lilly Date: Tue, 31 Oct 2023 15:53:35 -0700 Subject: [PATCH 17/34] Implement the second stage for the fine subcycling --- .../mpas_ocn_time_integration_fblts.F | 141 +++++++++++++++++- 1 file changed, 138 insertions(+), 3 deletions(-) diff --git a/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_fblts.F b/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_fblts.F index 6b0807e3d08e..8e5444f08c79 100644 --- a/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_fblts.F +++ b/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_fblts.F @@ -806,7 +806,7 @@ subroutine ocn_time_integrator_fblts(domain, dt)!{{{ ! Calculate predicted 'current' data at intermediate ! time-level on interface 1: - ! (u,h) = (im/M)*new + (1- (im/M))*current + ! u,h = (im/M)*new + (1- (im/M))*current call mpas_timer_start('FB_LTS interface prediction') normalVelocityForTend(:,:) = normalVelocityCur(:,:) @@ -952,7 +952,75 @@ subroutine ocn_time_integrator_fblts(domain, dt)!{{{ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! BEGIN: Thickness stage 2 ! Compute the second stage of the thickness on fine - + + ! Calculate predicted stage 1 data at intermediate + ! time-level on interface 1: + ! u,h = (im/M)*new + (1/M)*firstStage + (1-(im+1)/M)*current + call mpas_timer_start('FB_LTS interface prediction') + + normalVelocityForTend(:,:) = normalVelocityFirstStage(:,:) + do ie = 1, nEdgesInLTSRegion(1,2) + iEdge = edgesInLTSRegion(1,2,ie) + normalVelocityForTend(:,iEdge) = (im/M) * normalVelocityNew(:,iEdge) & + + (1/M) * normalvelocityFirstStage(:,iEdge) & + + (1 - ((im+1)/M)) * normalVelocityCur(:,iEdge) + end do + + layerThicknessForTend(:,:) = layerThicknessFirstStage(:,:) + do ic = 1, nCellsInLTSRegion(1,2) + iCell = cellsInLTSRegion(1,2,ic) + layerThicknessForTend(:,iCell) = (im/M) * layerThicknessNew(:,iCell) & + + (1/M) * layerThicknessFirstStage(:,iCell) & + + (1 - ((im+1)/M)) * layerThicknessCur(:,iCell) + end do + + call mpas_timer_stop('FB_LTS interface prediction') + + ! Halo update before tend calculation + call mpas_timer_start("FB_LTS prognostic halo update") + call mpas_timer_start("FB_LTS thickness halo update") + + call mpas_dmpar_field_halo_exch(domain, 'normalVelocity', timeLevel=5) + call mpas_dmpar_field_halo_exch(domain, 'layerThickness', timeLevel=5) + + call mpas_timer_stop("FB_LTS thickness halo update") + call mpas_timer_stop("FB_LTS prognostic halo update") + + ! Compute fast tendencies for thickness + call mpas_timer_start("FB_LTS compute fast tendencies") + call mpas_timer_start("FB_LTS compute fast thickness tendencies") + + call ocn_lts_thick_tend(LTSPool, tendPool, & + normalVelocityForTend, layerThicknessForTend, & + 1, 1, 0, 0, 0) + + call mpas_timer_stop("FB_LTS compute fast thickness tendencies") + call mpas_timer_stop("FB_LTS compute fast tendencies") + + ! Advance thickness solution with first stage + call mpas_timer_start("FB_LTS advance solution") + call mpas_timer_start("FB_LTS advance thickness solution") + + ! Fine cells adjacent to interface 1 + do ic = 1, nCellsInLTSRegion(1,3) + iCell = cellsInLTSRegion(1,3,ic) + do k = 1, maxLevelCell(iCell) + layerThicknessSecondStage(k,iCell) = layerThicknessCur(k,iCell) & + + weight2nd * dtFine * layerThicknessTend(k,iCell) + end do + end do + ! Fine in the interior, away from interface 1 + do ic = 1, nCellsInLTSRegion(1,1) + iCell = cellsInLTSRegion(1,1,ic) + do k = 1, maxLevelCell(iCell) + layerThicknessSecondStage(k,iCell) = layerThicknessCur(k,iCell) & + + weight2nd * dtFine * layerThicknessTend(k,iCell) + end do + end do + + call mpas_timer_stop("FB_LTS advance thickness solution") + call mpas_timer_stop("FB_LTS advance solution") + ! END: Thickness stage 2 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! @@ -960,7 +1028,74 @@ subroutine ocn_time_integrator_fblts(domain, dt)!{{{ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! BEGIN: Normal velocity stage 2 ! Compute the second stage of the normal velocity on fine - + + ! Perform forward-backward average of thickness for stage two of + ! FB-RK(3,2): + ! h = weight2*stage2 + (1-weight2)*current + layerThicknessForTend(:,:) = config_fb_weight_2 * layerThicknessSecondStage(:,:) & + + (1 - config_fb_weight_2) * layerThicknessCur(:,:) + + ! Calculate predicted FB averaged thickness data at + ! intermediate time-level on interface 1: + ! h = weight2*'secondStage' + (1-weight2)*'current' + ! = weight2*( (im/M)*new + (1/M)*secondStage + (1-(im+1)/M)*current ) + ! + (1-weight2)*( (im/M)*new + (1-im/M)*current ) + call mpas_timer_start('FB_LTS interface prediction') + + do ic = 1, nCellsInLTSRegion(1,2) + iCell = cellsInLTSRegion(1,2,ic) + layerThicknessForTend(:,iCell) = config_fb_weight_2 & + * ( (im/M) * layerThicknessNew(:,iCell) & + + (1/M) * layerThicknessSecondStage(:,iCell) & + + (1 - ((im+1)/M)) * layerThicknessCur(:,iCell) ) & + + (1 - config_fb_weight_2) & + * ( (im/M) * layerThicknessNew(:,iCell) & + + (1 - (im/M)) * layerThicknessCur(:,iCell) ) + end do + + call mpas_timer_stop('FB_LTS interface prediction') + + ! Halo update before tend calculation + call mpas_timer_start("FB_LTS prognostic halo update") + call mpas_timer_start("FB_LTS thickness halo update") + + ! Don't need to update normalVelocityForTend, previously updated + call mpas_dmpar_field_halo_exch(domain, 'layerThickness', timeLevel=5) + + call mpas_timer_stop("FB_LTS thickness halo update") + call mpas_timer_stop("FB_LTS prognostic halo update") + + ! Compute fast tendencies for normal velocity + call mpas_timer_start("FB_LTS compute fast tendencies") + call mpas_timer_start("FB_LTS compute fast velocity tendencies") + + call ocn_lts_vel_tend(LTSPool, tendPool, & + normalVelocityForTend, layerThicknessForTend, & + 1, 1, 0, 0, 0) + + call mpas_timer_stop("FB_LTS compute fast velocity tendencies") + call mpas_timer_stop("FB_LTS compute fast tendencies") + + ! Advance normal velocity solution with second stage + call mpas_timer_start("FB_LTS advance solution") + call mpas_timer_start("FB_LTS advance velocity solution") + + ! Fine cells adjacent to interface 1 + do ie = 1, nEdgesInLTSRegion(1,3) + iEdge = edgesInLTSRegion(1,3,ie) + normalVelocitySecondStage(:,iEdge) = normalVelocityCur(:,iEdge) & + + weight2nd * dtFine * normalVelocityTend(:,iEdge) + end do + ! Fine in the interior, away from interface 1 + do ie = 1, nEdgesInLTSRegion(1,1) + iEdge = edgesInLTSRegion(1,1,ie) + normalVelocitySecondStage(:,iEdge) = normalVelocityCur(:,iEdge) & + + weight2nd * dtFine * normalVelocityTend(:,iEdge) + end do + + call mpas_timer_stop("FB_LTS advance velocity solution") + call mpas_timer_stop("FB_LTS advance solution") + ! END: Normal velocity stage 2 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! From a0c6d253dda94f7b7c39722ffeaf53d268a78d35 Mon Sep 17 00:00:00 2001 From: Jeremy Lilly Date: Tue, 31 Oct 2023 16:39:21 -0700 Subject: [PATCH 18/34] Implement the third stage for the fine subcycling, also make loops over nVertLevels implicit for layerThickness advancements to match normalVelocity advancements --- .../mpas_ocn_time_integration_fblts.F | 240 +++++++++++++----- 1 file changed, 183 insertions(+), 57 deletions(-) diff --git a/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_fblts.F b/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_fblts.F index 8e5444f08c79..0968b26e1136 100644 --- a/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_fblts.F +++ b/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_fblts.F @@ -424,27 +424,21 @@ subroutine ocn_time_integrator_fblts(domain, dt)!{{{ do iRegion = 1,nRegions do ic = 1, nCellsInLTSRegion(iRegion,2) iCell = cellsInLTSRegion(iRegion,2,ic) - do k = 1, maxLevelCell(iCell) - layerThicknessFirstStage(k,iCell) = layerThicknessCur(k,iCell) & - + weight1st * dt * layerThicknessTend(k,iCell) - end do + layerThicknessFirstStage(:,iCell) = layerThicknessCur(:,iCell) & + + weight1st * dt * layerThicknessTend(:,iCell) end do end do ! Coarse do ic = 1, nCellsInLTSRegion(2,1) iCell = cellsInLTSRegion(2,1,ic) - do k = 1, maxLevelCell(iCell) - layerThicknessFirstStage(k,iCell) = layerThicknessCur(k,iCell) & - + weight1st * dt * layerThicknessTend(k,iCell) - end do + layerThicknessFirstStage(:,iCell) = layerThicknessCur(:,iCell) & + + weight1st * dt * layerThicknessTend(:,iCell) end do ! Fine layers close to interface layers do ic = 1, nCellsInLTSRegion(1,3) iCell = cellsInLTSRegion(1,3,ic) - do k = 1, maxLevelCell(iCell) - layerThicknessFirstStage(k,iCell) = layerThicknessCur(k,iCell) & - + weight1st * dt * layerThicknessTend(k,iCell) - end do + layerThicknessFirstStage(:,iCell) = layerThicknessCur(:,iCell) & + + weight1st * dt * layerThicknessTend(:,iCell) end do call mpas_timer_stop("FB_LTS advance thickness solution") @@ -556,27 +550,21 @@ subroutine ocn_time_integrator_fblts(domain, dt)!{{{ do iRegion = 1,nRegions do ic = 1, nCellsInLTSRegion(iRegion,2) iCell = cellsInLTSRegion(iRegion,2,ic) - do k = 1, maxLevelCell(iCell) - layerThicknessSecondStage(k,iCell) = layerThicknessCur(k,iCell) & - + weight2nd * dt * layerThicknessTend(k,iCell) - end do + layerThicknessSecondStage(:,iCell) = layerThicknessCur(:,iCell) & + + weight2nd * dt * layerThicknessTend(:,iCell) end do end do ! Coarse do ic = 1, nCellsInLTSRegion(2,1) iCell = cellsInLTSRegion(2,1,ic) - do k = 1, maxLevelCell(iCell) - layerThicknessSecondStage(k,iCell) = layerThicknessCur(k,iCell) & - + weight2nd * dt * layerThicknessTend(k,iCell) - end do + layerThicknessSecondStage(:,iCell) = layerThicknessCur(:,iCell) & + + weight2nd * dt * layerThicknessTend(:,iCell) end do ! Fine layers close to interface layers do ic = 1, nCellsInLTSRegion(1,3) iCell = cellsInLTSRegion(1,3,ic) - do k = 1, maxLevelCell(iCell) - layerThicknessSecondStage(k,iCell) = layerThicknessCur(k,iCell) & - + weight2nd * dt * layerThicknessTend(k,iCell) - end do + layerThicknessSecondStage(:,iCell) = layerThicknessCur(:,iCell) & + + weight2nd * dt * layerThicknessTend(:,iCell) end do call mpas_timer_stop("FB_LTS advance thickness solution") @@ -688,27 +676,21 @@ subroutine ocn_time_integrator_fblts(domain, dt)!{{{ do iRegion = 1,nRegions do ic = 1, nCellsInLTSRegion(iRegion,2) iCell = cellsInLTSRegion(iRegion,2,ic) - do k = 1, maxLevelCell(iCell) - layerThicknessNew(k,iCell) = layerThicknessCur(k,iCell) & - + dt * layerThicknessTend(k,iCell) - end do + layerThicknessNew(:,iCell) = layerThicknessCur(:,iCell) & + + dt * layerThicknessTend(:,iCell) end do end do ! Coarse do ic = 1, nCellsInLTSRegion(2,1) iCell = cellsInLTSRegion(2,1,ic) - do k = 1, maxLevelCell(iCell) - layerThicknessNew(k,iCell) = layerThicknessCur(k,iCell) & - + dt * layerThicknessTend(k,iCell) - end do + layerThicknessNew(:,iCell) = layerThicknessCur(:,iCell) & + + dt * layerThicknessTend(:,iCell) end do ! Fine layers close to interface layers do ic = 1, nCellsInLTSRegion(1,3) iCell = cellsInLTSRegion(1,3,ic) - do k = 1, maxLevelCell(iCell) - layerThicknessNew(k,iCell) = layerThicknessCur(k,iCell) & - + dt * layerThicknessTend(k,iCell) - end do + layerThicknessNew(:,iCell) = layerThicknessCur(:,iCell) & + + dt * layerThicknessTend(:,iCell) end do call mpas_timer_stop("FB_LTS advance thickness solution") @@ -853,18 +835,14 @@ subroutine ocn_time_integrator_fblts(domain, dt)!{{{ ! Fine cells adjacent to interface 1 do ic = 1, nCellsInLTSRegion(1,3) iCell = cellsInLTSRegion(1,3,ic) - do k = 1, maxLevelCell(iCell) - layerThicknessFirstStage(k,iCell) = layerThicknessCur(k,iCell) & - + weight1st * dtFine * layerThicknessTend(k,iCell) - end do + layerThicknessFirstStage(:,iCell) = layerThicknessCur(:,iCell) & + + weight1st * dtFine * layerThicknessTend(:,iCell) end do ! Fine in the interior, away from interface 1 do ic = 1, nCellsInLTSRegion(1,1) iCell = cellsInLTSRegion(1,1,ic) - do k = 1, maxLevelCell(iCell) - layerThicknessFirstStage(k,iCell) = layerThicknessCur(k,iCell) & - + weight1st * dtFine * layerThicknessTend(k,iCell) - end do + layerThicknessFirstStage(:,iCell) = layerThicknessCur(:,iCell) & + + weight1st * dtFine * layerThicknessTend(:,iCell) end do call mpas_timer_stop("FB_LTS advance thickness solution") @@ -962,7 +940,7 @@ subroutine ocn_time_integrator_fblts(domain, dt)!{{{ do ie = 1, nEdgesInLTSRegion(1,2) iEdge = edgesInLTSRegion(1,2,ie) normalVelocityForTend(:,iEdge) = (im/M) * normalVelocityNew(:,iEdge) & - + (1/M) * normalvelocityFirstStage(:,iEdge) & + + (1/M) * normalVelocityFirstStage(:,iEdge) & + (1 - ((im+1)/M)) * normalVelocityCur(:,iEdge) end do @@ -997,25 +975,21 @@ subroutine ocn_time_integrator_fblts(domain, dt)!{{{ call mpas_timer_stop("FB_LTS compute fast thickness tendencies") call mpas_timer_stop("FB_LTS compute fast tendencies") - ! Advance thickness solution with first stage + ! Advance thickness solution with second stage call mpas_timer_start("FB_LTS advance solution") call mpas_timer_start("FB_LTS advance thickness solution") ! Fine cells adjacent to interface 1 do ic = 1, nCellsInLTSRegion(1,3) iCell = cellsInLTSRegion(1,3,ic) - do k = 1, maxLevelCell(iCell) - layerThicknessSecondStage(k,iCell) = layerThicknessCur(k,iCell) & - + weight2nd * dtFine * layerThicknessTend(k,iCell) - end do + layerThicknessSecondStage(:,iCell) = layerThicknessCur(:,iCell) & + + weight2nd * dtFine * layerThicknessTend(:,iCell) end do ! Fine in the interior, away from interface 1 do ic = 1, nCellsInLTSRegion(1,1) iCell = cellsInLTSRegion(1,1,ic) - do k = 1, maxLevelCell(iCell) - layerThicknessSecondStage(k,iCell) = layerThicknessCur(k,iCell) & - + weight2nd * dtFine * layerThicknessTend(k,iCell) - end do + layerThicknessSecondStage(:,iCell) = layerThicknessCur(:,iCell) & + + weight2nd * dtFine * layerThicknessTend(:,iCell) end do call mpas_timer_stop("FB_LTS advance thickness solution") @@ -1104,7 +1078,79 @@ subroutine ocn_time_integrator_fblts(domain, dt)!{{{ ! BEGIN: Thickness stage 3 ! Compute the third stage of the thickness on fine and ! increment the third stage interface correction terms - + + ! Calculate predicted stage 2 data at intermediate + ! time-level on interface 1: + ! u,h = (im/M)*new + (1/M)*secondStage + (1-(im+1)/M)*current + call mpas_timer_start('FB_LTS interface prediction') + + normalVelocityForTend(:,:) = normalVelocitySecondStage(:,:) + do ie = 1, nEdgesInLTSRegion(1,2) + iEdge = edgesInLTSRegion(1,2,ie) + normalVelocityForTend(:,iEdge) = (im/M) * normalVelocityNew(:,iEdge) & + + (1/M) * normalVelocitySecondStage(:,iEdge) & + + (1 - ((im+1)/M)) * normalVelocityCur(:,iEdge) + end do + + layerThicknessForTend(:,:) = layerThicknessSecondStage(:,:) + do ic = 1, nCellsInLTSRegion(1,2) + iCell = cellsInLTSRegion(1,2,ic) + layerThicknessForTend(:,iCell) = (im/M) * layerThicknessNew(:,iCell) & + + (1/M) * layerThicknessSecondStage(:,iCell) & + + (1 - ((im+1)/M)) * layerThicknessCur(:,iCell) + end do + + call mpas_timer_stop('FB_LTS interface prediction') + + ! Halo update before tend calculation + call mpas_timer_start("FB_LTS prognostic halo update") + call mpas_timer_start("FB_LTS thickness halo update") + + call mpas_dmpar_field_halo_exch(domain, 'normalVelocity', timeLevel=5) + call mpas_dmpar_field_halo_exch(domain, 'layerThickness', timeLevel=5) + + call mpas_timer_stop("FB_LTS thickness halo update") + call mpas_timer_stop("FB_LTS prognostic halo update") + + ! Compute fast tendencies for thickness + call mpas_timer_start("FB_LTS compute fast tendencies") + call mpas_timer_start("FB_LTS compute fast thickness tendencies") + + call ocn_lts_thick_tend(LTSPool, tendPool, & + normalVelocityForTend, layerThicknessForTend, & + 1, 1, 1, 1, 0) + + call mpas_timer_stop("FB_LTS compute fast thickness tendencies") + call mpas_timer_stop("FB_LTS compute fast tendencies") + + ! Advance thickness solution with third stage + call mpas_timer_start("FB_LTS advance solution") + call mpas_timer_start("FB_LTS advance thickness solution") + + ! Increment interface correction terms + do iRegion = 1,nRegions + do ic = 1, nCellsInLTSRegion(iRegion,2) + iCell = cellsInLTSRegion(iRegion,2,ic) + layerThicknessTendSum3rd(:,iCell) = layerThicknessTendSum3rd(:,iCell) & + + layerThicknessTend(:,iCell) + end do + end do + ! Fine cells adjacent to interface 1 + do ic = 1, nCellsInLTSRegion(1,3) + iCell = cellsInLTSRegion(1,3,ic) + layerThicknessNew(:,iCell) = layerThicknessCur(:,iCell) & + + dtFine * layerThicknessTend(:,iCell) + end do + ! Fine in the interior, away from interface 1 + do ic = 1, nCellsInLTSRegion(1,1) + iCell = cellsInLTSRegion(1,1,ic) + layerThicknessNew(:,iCell) = layerThicknessCur(:,iCell) & + + dtFine * layerThicknessTend(:,iCell) + end do + + call mpas_timer_stop("FB_LTS advance thickness solution") + call mpas_timer_stop("FB_LTS advance solution") + ! END: Thickness stage 3 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! @@ -1112,8 +1158,88 @@ subroutine ocn_time_integrator_fblts(domain, dt)!{{{ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! BEGIN: Normal velocity stage 3 ! Compute the third stage of the normal velocity on fine and - ! increment the thrid stage interface correction terms - + ! increment the third stage interface correction terms + + ! Perform forward-backward average of thickness for stage three of + ! FB-RK(3,2): + ! h = weight3*new + (1-2*weight3)*stageTwo + weight3*current + layerThicknessForTend(:,:) = config_fb_weight_3 * layerThicknessNew(:,:) & + + (1 - 2*config_fb_weight_3) * layerThicknessSecondStage(:,:) & + + config_fb_weight_3 * layerThicknessCur(:,:) + + ! Calculate predicted FB averaged thickness data at + ! intermediate time-level on interface 1: + ! h = weight3*new + (1-2*weight3)*'secondStage' + weight3*'current' + ! = weight3*( ((im+1)/M)*new + (1-(im+1)/M)*current ) + ! + (1-2*weight3)*( (im/M)*new + (1/M)*secondStage + (1-(im+1)/M)*current ) + ! + weight3*( (im/M)*new + (1-im/M)*current ) + call mpas_timer_start('FB_LTS interface prediction') + + do ic = 1, nCellsInLTSRegion(1,2) + iCell = cellsInLTSRegion(1,2,ic) + layerThicknessForTend(:,iCell) = config_fb_weight_3 & + * ( ((im+1)/M) * layerThicknessNew(:,iCell) & + + (1 - ((im+1)/M)) * layerThicknessCur(:,iCell) ) & + + (1 - 2*config_fb_weight_3) & + * ( (im/M) * layerThicknessNew(:,iCell) & + + (1/M) * layerThicknessSecondStage(:,iCell) & + + (1 - ((im+1)/M)) * layerThicknessCur(:,iCell) ) & + + config_fb_weight_3 & + * ( (im/M) * layerThicknessNew(:,iCell) & + + (1 - (im/M)) * layerThicknessCur(:,iCell) ) + end do + + call mpas_timer_stop('FB_LTS interface prediction') + + ! Halo update before tend calculation + call mpas_timer_start("FB_LTS prognostic halo update") + call mpas_timer_start("FB_LTS thickness halo update") + + ! Don't need to update normalVelocityForTend, previously updated + call mpas_dmpar_field_halo_exch(domain, 'layerThickness', timeLevel=5) + + call mpas_timer_stop("FB_LTS thickness halo update") + call mpas_timer_stop("FB_LTS prognostic halo update") + + ! Compute fast tendencies for normal velocity + call mpas_timer_start("FB_LTS compute fast tendencies") + call mpas_timer_start("FB_LTS compute fast velocity tendencies") + + call ocn_lts_vel_tend(LTSPool, tendPool, & + normalVelocityForTend, layerThicknessForTend, & + 1, 1, 1, 1, 0) + + call mpas_timer_stop("FB_LTS compute fast velocity tendencies") + call mpas_timer_stop("FB_LTS compute fast tendencies") + + ! Advance normal velocity solution with third stage + call mpas_timer_start("FB_LTS advance solution") + call mpas_timer_start("FB_LTS advance velocity solution") + + ! Increment interface correction terms + do iRegion = 1,nRegions + do ie = 1, nEdgesInLTSRegion(iRegion,2) + iEdge = edgesInLTSRegion(iRegion,2,ie) + normalVelocityTendSum3rd(:,iEdge) = normalVelocityTendSum3rd(:,iEdge) & + + normalVelocityTend(:,iEdge) + end do + end do + ! Fine cells adjacent to interface 1 + do ie = 1, nEdgesInLTSRegion(1,3) + iEdge = edgesInLTSRegion(1,3,ie) + normalVelocityNew(:,iEdge) = normalVelocityCur(:,iEdge) & + + dtFine * normalVelocityTend(:,iEdge) + end do + ! Fine in the interior, away from interface 1 + do ie = 1, nEdgesInLTSRegion(1,1) + iEdge = edgesInLTSRegion(1,1,ie) + normalVelocityNew(:,iEdge) = normalVelocityCur(:,iEdge) & + + dtFine * normalVelocityTend(:,iEdge) + end do + + call mpas_timer_stop("FB_LTS advance velocity solution") + call mpas_timer_stop("FB_LTS advance solution") + ! END: Normal velocity stage 3 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! From e54084b8248fe7de631e1863478473e7a8e3dd43 Mon Sep 17 00:00:00 2001 From: Jeremy Lilly Date: Tue, 31 Oct 2023 16:52:27 -0700 Subject: [PATCH 19/34] Implement the interface correction step --- .../mpas_ocn_time_integration_fblts.F | 33 +++++++++++++++++++ 1 file changed, 33 insertions(+) diff --git a/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_fblts.F b/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_fblts.F index 0968b26e1136..836cf896ad90 100644 --- a/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_fblts.F +++ b/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_fblts.F @@ -1253,10 +1253,43 @@ subroutine ocn_time_integrator_fblts(domain, dt)!{{{ ! BEGIN INTERFACE CORRECTION !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! Calculate the corrected solution in the interface layers + call mpas_timer_start("FB_LTS interface correction") + + ! Perform the correction on interface 1 and 2 + do iRegion = 1, nRegions + ! Normal velocity correction + do ie = 1, nEdgesInLTSRegion(iRegion,2) + iEdge = edgesInLTSRegion(iRegion,2,ie) + normalVelocityNew(:,iEdge) = normalVelocityCur(:,iEdge) & + + dtFine * normalVelocityTendSum3rd(:,iEdge) + end do + ! Layer thickness correction + do ic = 1, nCellsInLTSRegion(iRegion,2) + iCell = cellsInLTSRegion(iRegion,2,ic) + layerThicknessNew(:,iCell) = layerThicknessCur(:,iCell) & + + dtFine * layerThicknessTendSum3rd(:,iCell) + end do + end do + + call mpas_timer_stop("FB_LTS interface correction") !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! END INTERFACE CORRECTION !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + ! Final prognostic halo update + call mpas_timer_start("FB_LTS prognostic halo update") + + call mpas_timer_start("FB_LTS velocity halo update") + call mpas_dmpar_field_halo_exch(domain, 'normalVelocity', timeLevel=2) + call mpas_timer_stop("FB_LTS velocity halo update") + + call mpas_timer_start("FB_LTS thickness halo update") + call mpas_dmpar_field_halo_exch(domain, 'layerThickness', timeLevel=2) + call mpas_timer_stop("FB_LTS thickness halo update") + + call mpas_timer_stop("FB_LTS prognostic halo update") call mpas_timer_stop("FB_LTS main loop") From 81e273f2f12f87343b71178860bd5d8babcba763 Mon Sep 17 00:00:00 2001 From: Jeremy Lilly Date: Tue, 31 Oct 2023 17:01:16 -0700 Subject: [PATCH 20/34] Move timers around to correctly differentiate between halo updates for thickness vs velocity --- .../mpas_ocn_time_integration_fblts.F | 89 ++++++++++++------- 1 file changed, 55 insertions(+), 34 deletions(-) diff --git a/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_fblts.F b/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_fblts.F index 836cf896ad90..c753a4492303 100644 --- a/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_fblts.F +++ b/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_fblts.F @@ -462,12 +462,13 @@ subroutine ocn_time_integrator_fblts(domain, dt)!{{{ ! Halo update before tend calculation ! Might be able to remove this? call mpas_timer_start("FB_LTS prognostic halo update") - call mpas_timer_start("FB_LTS thickness halo update") - + ! Don't need to update normalVelocityCur, already updated + + call mpas_timer_start("FB_LTS thickness halo update") call mpas_dmpar_field_halo_exch(domain, 'layerThickness', timeLevel=5) - call mpas_timer_stop("FB_LTS thickness halo update") + call mpas_timer_stop("FB_LTS prognostic halo update") ! Compute fast tendencies for normal velocity @@ -523,12 +524,15 @@ subroutine ocn_time_integrator_fblts(domain, dt)!{{{ ! Halo update before tend calculation call mpas_timer_start("FB_LTS prognostic halo update") + call mpas_timer_start("FB_LTS velocity halo update") - call mpas_dmpar_field_halo_exch(domain, 'normalVelocity', timeLevel=3) - call mpas_dmpar_field_halo_exch(domain, 'layerThickness', timeLevel=3) - call mpas_timer_stop("FB_LTS velocity halo update") + + call mpas_timer_start("FB_LTS thickness halo update") + call mpas_dmpar_field_halo_exch(domain, 'layerThickness', timeLevel=3) + call mpas_timer_stop("FB_LTS thickness halo update") + call mpas_timer_stop("FB_LTS prognostic halo update") ! Compute fast tendencies for thickness @@ -588,12 +592,13 @@ subroutine ocn_time_integrator_fblts(domain, dt)!{{{ ! Halo update before tend calculation ! Might be able to remove this? call mpas_timer_start("FB_LTS prognostic halo update") - call mpas_timer_start("FB_LTS thickness halo update") - + ! Don't need to update normalVelocityFirstStage, previously updated + + call mpas_timer_start("FB_LTS thickness halo update") call mpas_dmpar_field_halo_exch(domain, 'layerThickness', timeLevel=5) - call mpas_timer_stop("FB_LTS thickness halo update") + call mpas_timer_stop("FB_LTS prognostic halo update") ! Compute fast tendencies for normal velocity @@ -649,12 +654,15 @@ subroutine ocn_time_integrator_fblts(domain, dt)!{{{ ! Halo update before tends calculation call mpas_timer_start("FB_LTS prognostic halo update") + call mpas_timer_start("FB_LTS velocity halo update") - call mpas_dmpar_field_halo_exch(domain, 'normalVelocity', timeLevel=4) - call mpas_dmpar_field_halo_exch(domain, 'layerThickness', timeLevel=4) - call mpas_timer_stop("FB_LTS velocity halo update") + + call mpas_timer_start("FB_LTS thickness halo update") + call mpas_dmpar_field_halo_exch(domain, 'layerThickness', timeLevel=4) + call mpas_timer_stop("FB_LTS thickness halo update") + call mpas_timer_stop("FB_LTS prognostic halo update") ! Compute fast tendencies for thickness @@ -715,12 +723,13 @@ subroutine ocn_time_integrator_fblts(domain, dt)!{{{ ! Halo update before tend calculation ! Might be able to remove this? call mpas_timer_start("FB_LTS prognostic halo update") - call mpas_timer_start("FB_LTS thickness halo update") - + ! Don't need to update normalVelocitySecondStage, previously updated + + call mpas_timer_start("FB_LTS thickness halo update") call mpas_dmpar_field_halo_exch(domain, 'layerThickness', timeLevel=5) - call mpas_timer_stop("FB_LTS thickness halo update") + call mpas_timer_stop("FB_LTS prognostic halo update") ! Compute fast tendencies for normal velocity @@ -762,11 +771,11 @@ subroutine ocn_time_integrator_fblts(domain, dt)!{{{ ! Halo update ! Might be able to remove this? call mpas_timer_start("FB_LTS prognostic halo update") + call mpas_timer_start("FB_LTS velocity halo update") - call mpas_dmpar_field_halo_exch(domain, 'normalVelocity', timeLevel=2) - call mpas_timer_stop("FB_LTS velocity halo update") + call mpas_timer_stop("FB_LTS prognostic halo update") ! END: Normal velocity stage 3 @@ -809,12 +818,15 @@ subroutine ocn_time_integrator_fblts(domain, dt)!{{{ ! Halo update before tend calculation call mpas_timer_start("FB_LTS prognostic halo update") - call mpas_timer_start("FB_LTS thickness halo update") - + + call mpas_timer_start("FB_LTS velocity halo update") call mpas_dmpar_field_halo_exch(domain, 'normalVelocity', timeLevel=5) + call mpas_timer_stop("FB_LTS velocity halo update") + + call mpas_timer_start("FB_LTS thickness halo update") call mpas_dmpar_field_halo_exch(domain, 'layerThickness', timeLevel=5) - call mpas_timer_stop("FB_LTS thickness halo update") + call mpas_timer_stop("FB_LTS prognostic halo update") ! Compute fast tendencies for thickness @@ -884,12 +896,13 @@ subroutine ocn_time_integrator_fblts(domain, dt)!{{{ ! Halo update before tend calculation call mpas_timer_start("FB_LTS prognostic halo update") - call mpas_timer_start("FB_LTS thickness halo update") - + ! Don't need to update normalVelocityForTend, previously updated + + call mpas_timer_start("FB_LTS thickness halo update") call mpas_dmpar_field_halo_exch(domain, 'layerThickness', timeLevel=5) - call mpas_timer_stop("FB_LTS thickness halo update") + call mpas_timer_stop("FB_LTS prognostic halo update") ! Compute fast tendencies for normal velocity @@ -956,12 +969,15 @@ subroutine ocn_time_integrator_fblts(domain, dt)!{{{ ! Halo update before tend calculation call mpas_timer_start("FB_LTS prognostic halo update") - call mpas_timer_start("FB_LTS thickness halo update") - + + call mpas_timer_start("FB_LTS velocity halo update") call mpas_dmpar_field_halo_exch(domain, 'normalVelocity', timeLevel=5) + call mpas_timer_stop("FB_LTS velocity halo update") + + call mpas_timer_start("FB_LTS thickness halo update") call mpas_dmpar_field_halo_exch(domain, 'layerThickness', timeLevel=5) - call mpas_timer_stop("FB_LTS thickness halo update") + call mpas_timer_stop("FB_LTS prognostic halo update") ! Compute fast tendencies for thickness @@ -1031,12 +1047,13 @@ subroutine ocn_time_integrator_fblts(domain, dt)!{{{ ! Halo update before tend calculation call mpas_timer_start("FB_LTS prognostic halo update") - call mpas_timer_start("FB_LTS thickness halo update") - + ! Don't need to update normalVelocityForTend, previously updated + + call mpas_timer_start("FB_LTS thickness halo update") call mpas_dmpar_field_halo_exch(domain, 'layerThickness', timeLevel=5) - call mpas_timer_stop("FB_LTS thickness halo update") + call mpas_timer_stop("FB_LTS prognostic halo update") ! Compute fast tendencies for normal velocity @@ -1104,12 +1121,15 @@ subroutine ocn_time_integrator_fblts(domain, dt)!{{{ ! Halo update before tend calculation call mpas_timer_start("FB_LTS prognostic halo update") - call mpas_timer_start("FB_LTS thickness halo update") - + + call mpas_timer_start("FB_LTS velocity halo update") call mpas_dmpar_field_halo_exch(domain, 'normalVelocity', timeLevel=5) + call mpas_timer_stop("FB_LTS velocity halo update") + + call mpas_timer_start("FB_LTS thickness halo update") call mpas_dmpar_field_halo_exch(domain, 'layerThickness', timeLevel=5) - call mpas_timer_stop("FB_LTS thickness halo update") + call mpas_timer_stop("FB_LTS prognostic halo update") ! Compute fast tendencies for thickness @@ -1193,12 +1213,13 @@ subroutine ocn_time_integrator_fblts(domain, dt)!{{{ ! Halo update before tend calculation call mpas_timer_start("FB_LTS prognostic halo update") - call mpas_timer_start("FB_LTS thickness halo update") ! Don't need to update normalVelocityForTend, previously updated + + call mpas_timer_start("FB_LTS thickness halo update") call mpas_dmpar_field_halo_exch(domain, 'layerThickness', timeLevel=5) - call mpas_timer_stop("FB_LTS thickness halo update") + call mpas_timer_stop("FB_LTS prognostic halo update") ! Compute fast tendencies for normal velocity From c9b519671c5b10e165e314e6fe0962dd6f450fe5 Mon Sep 17 00:00:00 2001 From: Jeremy Lilly Date: Tue, 31 Oct 2023 17:18:39 -0700 Subject: [PATCH 21/34] Implement the cleanup phase at the end of the FB_LTS routine --- .../mpas_ocn_time_integration_fblts.F | 67 +++++++++++++++++-- 1 file changed, 63 insertions(+), 4 deletions(-) diff --git a/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_fblts.F b/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_fblts.F index c753a4492303..ff4d8e58ba10 100644 --- a/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_fblts.F +++ b/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_fblts.F @@ -32,7 +32,7 @@ module ocn_time_integration_fblts use ocn_constants use ocn_tendency - !use ocn_diagnostics + use ocn_diagnostics use ocn_mesh use ocn_vmix use ocn_config @@ -1321,9 +1321,68 @@ subroutine ocn_time_integrator_fblts(domain, dt)!{{{ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - call mpas_timer_start("FB_LTS cleanup phase") - - call mpas_timer_stop("FB_LTS cleanup phase") + call mpas_timer_start("FB_LTS cleanup") + + if (config_prescribe_velocity) then + do iEdge = 1, nEdgesAll + normalVelocityNew(:,iEdge) = normalVelocityCur(:,iEdge) + end do + end if + + if (config_prescribe_thickness) then + do iCell = 1, nCellsAll + layerThicknessNew(:,iCell) = layerThicknessCur(:,iCell) + end do + end if + + do iEdge = 1, nEdgesAll + normalTransportVelocity(:,iEdge) = normalVelocityNew(:,iEdge) + end do + + call mpas_reconstruct(meshPool, normalVelocityNew, & + velocityX, velocityY, velocityZ, & + velocityZonal, velocityMeridional, & + includeHalos = .true.) + + call mpas_reconstruct(meshPool, gradSSH, & + gradSSHX, gradSSHY, gradSSHZ, & + gradSSHZonal, gradSSHMeridional, & + includeHalos = .true.) + + do iCell = 1, nCellsAll + surfaceVelocity(indexSurfaceVelocityZonal, iCell) = velocityZonal(1, iCell) + surfaceVelocity(indexSurfaceVelocityMeridional, iCell) = velocityMeridional(1, iCell) + + SSHGradient(indexSSHGradientZonal, iCell) = gradSSHZonal(iCell) + SSHGradient(indexSSHGradientMeridional, iCell) = gradSSHMeridional(iCell) + end do + + call ocn_time_average_coupled_accumulate(statePool, forcingPool, 2) + + do iCell = 1, nCellsAll + surfaceVelocity(indexSurfaceVelocityZonal, iCell) = velocityZonal(1, iCell) + surfaceVelocity(indexSurfaceVelocityMeridional, iCell) = velocityMeridional(1, iCell) + + SSHGradient(indexSSHGradientZonal, iCell) = gradSSHZonal(iCell) + SSHGradient(indexSSHGradientMeridional, iCell) = gradSSHMeridional(iCell) + end do + + + ! Update diagnostics + call ocn_diagnostic_solve(dt, statePool, forcingPool, meshPool, & + verticalMeshPool, scratchPool, tracersPool, 2) + + !call mpas_pool_destroy_pool(tendSum1stPool) + !call mpas_pool_destroy_pool(tendSum2ndPool) + call mpas_pool_destroy_pool(tendSum3rdPool) + call mpas_pool_destroy_pool(tendSlowPool) + + !call mpas_pool_remove_subpool(block % structs, 'tend_sum_1st') + !call mpas_pool_remove_subpool(block % structs, 'tend_sum_2nd') + call mpas_pool_remove_subpool(block % structs, 'tend_sum_3rd') + call mpas_pool_remove_subpool(block % structs, 'tend_slow') + + call mpas_timer_stop("FB_LTS cleanup") end subroutine ocn_time_integrator_fblts!}}} From 0a477003d4ec60141bc0ff15d561a8b0f7398898 Mon Sep 17 00:00:00 2001 From: Jeremy Lilly Date: Wed, 8 Nov 2023 10:34:55 -0800 Subject: [PATCH 22/34] Add slow tends to the velocity tends during the substepping --- .../mode_forward/mpas_ocn_time_integration_fblts.F | 11 ++++++++++- 1 file changed, 10 insertions(+), 1 deletion(-) diff --git a/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_fblts.F b/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_fblts.F index ff4d8e58ba10..7a35531f5773 100644 --- a/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_fblts.F +++ b/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_fblts.F @@ -915,11 +915,14 @@ subroutine ocn_time_integrator_fblts(domain, dt)!{{{ call mpas_timer_stop("FB_LTS compute fast velocity tendencies") call mpas_timer_stop("FB_LTS compute fast tendencies") - + ! Advance normal velocity solution with first stage call mpas_timer_start("FB_LTS advance solution") call mpas_timer_start("FB_LTS advance velocity solution") + normalVelocityTend(:,:) = normalVelocityTend(:,:) & + + normalVelocityTendSlow(:,:) + ! Fine cells adjacent to interface 1 do ie = 1, nEdgesInLTSRegion(1,3) iEdge = edgesInLTSRegion(1,3,ie) @@ -1071,6 +1074,9 @@ subroutine ocn_time_integrator_fblts(domain, dt)!{{{ call mpas_timer_start("FB_LTS advance solution") call mpas_timer_start("FB_LTS advance velocity solution") + normalVelocityTend(:,:) = normalVelocityTend(:,:) & + + normalVelocityTendSlow(:,:) + ! Fine cells adjacent to interface 1 do ie = 1, nEdgesInLTSRegion(1,3) iEdge = edgesInLTSRegion(1,3,ie) @@ -1236,6 +1242,9 @@ subroutine ocn_time_integrator_fblts(domain, dt)!{{{ ! Advance normal velocity solution with third stage call mpas_timer_start("FB_LTS advance solution") call mpas_timer_start("FB_LTS advance velocity solution") + + normalVelocityTend(:,:) = normalVelocityTend(:,:) & + + normalVelocityTendSlow(:,:) ! Increment interface correction terms do iRegion = 1,nRegions From b39db88133cae4cbf4ca6a653b2555a51633663f Mon Sep 17 00:00:00 2001 From: Jeremy Lilly Date: Fri, 17 Nov 2023 16:14:23 -0800 Subject: [PATCH 23/34] Make sure that literals are reals, not integers, and cast some integers are reals to avoid integer division --- .../mpas_ocn_time_integration_fblts.F | 85 ++++++++++--------- 1 file changed, 43 insertions(+), 42 deletions(-) diff --git a/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_fblts.F b/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_fblts.F index 7a35531f5773..772caa6b19e7 100644 --- a/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_fblts.F +++ b/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_fblts.F @@ -457,7 +457,7 @@ subroutine ocn_time_integrator_fblts(domain, dt)!{{{ ! FB-RK(3,2): ! h = weight1*stage1 + (1-weight1)*current layerThicknessForTend(:,:) = config_fb_weight_1 * layerThicknessFirstStage(:,:) & - + (1 - config_fb_weight_1) * layerThicknessCur(:,:) + + (1.0_RKIND - config_fb_weight_1) * layerThicknessCur(:,:) ! Halo update before tend calculation ! Might be able to remove this? @@ -587,7 +587,7 @@ subroutine ocn_time_integrator_fblts(domain, dt)!{{{ ! FB-RK(3,2): ! h = weight2*stage2 + (1-weight2)*current layerThicknessForTend(:,:) = config_fb_weight_2 * layerThicknessSecondStage(:,:) & - + (1 - config_fb_weight_2) * layerThicknessCur(:,:) + + (1.0_RKIND - config_fb_weight_2) * layerThicknessCur(:,:) ! Halo update before tend calculation ! Might be able to remove this? @@ -717,7 +717,7 @@ subroutine ocn_time_integrator_fblts(domain, dt)!{{{ ! FB-RK(3,2): ! h = weight3*new + (1-2*weight3)*secondStage + weight3*current layerThicknessForTend(:,:) = config_fb_weight_3 * layerThicknessNew(:,:) & - + (1 - 2*config_fb_weight_3) * layerThicknessSecondStage(:,:) & + + (1.0_RKIND - 2.0_RKIND*config_fb_weight_3) * layerThicknessSecondStage(:,:) & + config_fb_weight_3 * layerThicknessCur(:,:) ! Halo update before tend calculation @@ -791,6 +791,7 @@ subroutine ocn_time_integrator_fblts(domain, dt)!{{{ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! do im = 0, M-1 + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! BEGIN: Thickness stage 1 ! Compute the first stage of the thickness on fine @@ -803,15 +804,15 @@ subroutine ocn_time_integrator_fblts(domain, dt)!{{{ normalVelocityForTend(:,:) = normalVelocityCur(:,:) do ie = 1, nEdgesInLTSRegion(1,2) iEdge = edgesInLTSRegion(1,2,ie) - normalVelocityForTend(:,iEdge) = (im/M) * normalVelocityNew(:,iEdge) & - + (1 - (im/M)) * normalVelocityCur(:,iEdge) + normalVelocityForTend(:,iEdge) = (REAL(im)/M) * normalVelocityNew(:,iEdge) & + + (1.0_RKIND - REAL(im)/M) * normalVelocityCur(:,iEdge) end do layerThicknessForTend(:,:) = layerThicknessCur(:,:) do ic = 1, nCellsInLTSRegion(1,2) iCell = cellsInLTSRegion(1,2,ic) - layerThicknessForTend(:,iCell) = (im/M) * layerThicknessNew(:,iCell) & - + (1 - (im/M)) * layerThicknessCur(:,iCell) + layerThicknessForTend(:,iCell) = (REAL(im)/M) * layerThicknessNew(:,iCell) & + + (1.0_RKIND - REAL(im)/M) * layerThicknessCur(:,iCell) end do call mpas_timer_stop('FB_LTS interface prediction') @@ -872,7 +873,7 @@ subroutine ocn_time_integrator_fblts(domain, dt)!{{{ ! FB-RK(3,2): ! h = weight1*stage1 + (1-weight1)*current layerThicknessForTend(:,:) = config_fb_weight_1 * layerThicknessFirstStage(:,:) & - + (1 - config_fb_weight_1) * layerThicknessCur(:,:) + + (1.0_RKIND - config_fb_weight_1) * layerThicknessCur(:,:) ! Calculate predicted FB averaged thickness data at ! intermediate time-level on interface 1: @@ -884,12 +885,12 @@ subroutine ocn_time_integrator_fblts(domain, dt)!{{{ do ic = 1, nCellsInLTSRegion(1,2) iCell = cellsInLTSRegion(1,2,ic) layerThicknessForTend(:,iCell) = config_fb_weight_1 & - * ( (im/M) * layerThicknessNew(:,iCell) & - + (1/M) * layerThicknessFirstStage(:,iCell) & - + (1 - ((im+1)/M)) * layerThicknessCur(:,iCell) ) & - + (1 - config_fb_weight_1) & - * ( (im/M) * layerThicknessNew(:,iCell) & - + (1 - (im/M)) * layerThicknessCur(:,iCell) ) + * ( (REAL(im)/M) * layerThicknessNew(:,iCell) & + + (1.0_RKIND/M) * layerThicknessFirstStage(:,iCell) & + + (1.0_RKIND - (REAL(im)+1.0_RKIND)/M) * layerThicknessCur(:,iCell) ) & + + (1.0_RKIND - config_fb_weight_1) & + * ( (REAL(im)/M) * layerThicknessNew(:,iCell) & + + (1.0_RKIND - REAL(im)/M) * layerThicknessCur(:,iCell) ) end do call mpas_timer_stop('FB_LTS interface prediction') @@ -955,17 +956,17 @@ subroutine ocn_time_integrator_fblts(domain, dt)!{{{ normalVelocityForTend(:,:) = normalVelocityFirstStage(:,:) do ie = 1, nEdgesInLTSRegion(1,2) iEdge = edgesInLTSRegion(1,2,ie) - normalVelocityForTend(:,iEdge) = (im/M) * normalVelocityNew(:,iEdge) & - + (1/M) * normalVelocityFirstStage(:,iEdge) & - + (1 - ((im+1)/M)) * normalVelocityCur(:,iEdge) + normalVelocityForTend(:,iEdge) = (REAL(im)/M) * normalVelocityNew(:,iEdge) & + + (1.0_RKIND/M) * normalVelocityFirstStage(:,iEdge) & + + (1.0_RKIND - (REAL(im)+1.0_RKIND)/M) * normalVelocityCur(:,iEdge) end do layerThicknessForTend(:,:) = layerThicknessFirstStage(:,:) do ic = 1, nCellsInLTSRegion(1,2) iCell = cellsInLTSRegion(1,2,ic) - layerThicknessForTend(:,iCell) = (im/M) * layerThicknessNew(:,iCell) & - + (1/M) * layerThicknessFirstStage(:,iCell) & - + (1 - ((im+1)/M)) * layerThicknessCur(:,iCell) + layerThicknessForTend(:,iCell) = (REAL(im)/M) * layerThicknessNew(:,iCell) & + + (1.0_RKIND/M) * layerThicknessFirstStage(:,iCell) & + + (1.0_RKIND - (REAL(im)+1.0_RKIND)/M) * layerThicknessCur(:,iCell) end do call mpas_timer_stop('FB_LTS interface prediction') @@ -1026,7 +1027,7 @@ subroutine ocn_time_integrator_fblts(domain, dt)!{{{ ! FB-RK(3,2): ! h = weight2*stage2 + (1-weight2)*current layerThicknessForTend(:,:) = config_fb_weight_2 * layerThicknessSecondStage(:,:) & - + (1 - config_fb_weight_2) * layerThicknessCur(:,:) + + (1.0_RKIND - config_fb_weight_2) * layerThicknessCur(:,:) ! Calculate predicted FB averaged thickness data at ! intermediate time-level on interface 1: @@ -1038,12 +1039,12 @@ subroutine ocn_time_integrator_fblts(domain, dt)!{{{ do ic = 1, nCellsInLTSRegion(1,2) iCell = cellsInLTSRegion(1,2,ic) layerThicknessForTend(:,iCell) = config_fb_weight_2 & - * ( (im/M) * layerThicknessNew(:,iCell) & - + (1/M) * layerThicknessSecondStage(:,iCell) & - + (1 - ((im+1)/M)) * layerThicknessCur(:,iCell) ) & - + (1 - config_fb_weight_2) & - * ( (im/M) * layerThicknessNew(:,iCell) & - + (1 - (im/M)) * layerThicknessCur(:,iCell) ) + * ( (REAL(im)/M) * layerThicknessNew(:,iCell) & + + (1.0_RKIND/M) * layerThicknessSecondStage(:,iCell) & + + (1.0_RKIND - (REAL(im)+1.0_RKIND)/M) * layerThicknessCur(:,iCell) ) & + + (1.0_RKIND - config_fb_weight_2) & + * ( (REAL(im)/M) * layerThicknessNew(:,iCell) & + + (1.0_RKIND - REAL(im)/M) * layerThicknessCur(:,iCell) ) end do call mpas_timer_stop('FB_LTS interface prediction') @@ -1110,17 +1111,17 @@ subroutine ocn_time_integrator_fblts(domain, dt)!{{{ normalVelocityForTend(:,:) = normalVelocitySecondStage(:,:) do ie = 1, nEdgesInLTSRegion(1,2) iEdge = edgesInLTSRegion(1,2,ie) - normalVelocityForTend(:,iEdge) = (im/M) * normalVelocityNew(:,iEdge) & - + (1/M) * normalVelocitySecondStage(:,iEdge) & - + (1 - ((im+1)/M)) * normalVelocityCur(:,iEdge) + normalVelocityForTend(:,iEdge) = (REAL(im)/M) * normalVelocityNew(:,iEdge) & + + (1.0_RKIND/M) * normalVelocitySecondStage(:,iEdge) & + + (1.0_RKIND - (REAL(im)+1.0_RKIND)/M) * normalVelocityCur(:,iEdge) end do layerThicknessForTend(:,:) = layerThicknessSecondStage(:,:) do ic = 1, nCellsInLTSRegion(1,2) iCell = cellsInLTSRegion(1,2,ic) - layerThicknessForTend(:,iCell) = (im/M) * layerThicknessNew(:,iCell) & - + (1/M) * layerThicknessSecondStage(:,iCell) & - + (1 - ((im+1)/M)) * layerThicknessCur(:,iCell) + layerThicknessForTend(:,iCell) = (REAL(im)/M) * layerThicknessNew(:,iCell) & + + (1.0_RKIND/M) * layerThicknessSecondStage(:,iCell) & + + (1.0_RKIND - (REAL(im)+1.0_RKIND)/M) * layerThicknessCur(:,iCell) end do call mpas_timer_stop('FB_LTS interface prediction') @@ -1190,7 +1191,7 @@ subroutine ocn_time_integrator_fblts(domain, dt)!{{{ ! FB-RK(3,2): ! h = weight3*new + (1-2*weight3)*stageTwo + weight3*current layerThicknessForTend(:,:) = config_fb_weight_3 * layerThicknessNew(:,:) & - + (1 - 2*config_fb_weight_3) * layerThicknessSecondStage(:,:) & + + (1.0_RKIND - 2.0_RKIND*config_fb_weight_3) * layerThicknessSecondStage(:,:) & + config_fb_weight_3 * layerThicknessCur(:,:) ! Calculate predicted FB averaged thickness data at @@ -1204,15 +1205,15 @@ subroutine ocn_time_integrator_fblts(domain, dt)!{{{ do ic = 1, nCellsInLTSRegion(1,2) iCell = cellsInLTSRegion(1,2,ic) layerThicknessForTend(:,iCell) = config_fb_weight_3 & - * ( ((im+1)/M) * layerThicknessNew(:,iCell) & - + (1 - ((im+1)/M)) * layerThicknessCur(:,iCell) ) & - + (1 - 2*config_fb_weight_3) & - * ( (im/M) * layerThicknessNew(:,iCell) & - + (1/M) * layerThicknessSecondStage(:,iCell) & - + (1 - ((im+1)/M)) * layerThicknessCur(:,iCell) ) & + * ( ((REAL(im)+1.0_RKIND)/M) * layerThicknessNew(:,iCell) & + + (1.0_RKIND - (REAL(im)+1.0_RKIND)/M) * layerThicknessCur(:,iCell) ) & + + (1.0_RKIND - 2.0_RKIND*config_fb_weight_3) & + * ( (REAL(im)/M) * layerThicknessNew(:,iCell) & + + (1.0_RKIND/M) * layerThicknessSecondStage(:,iCell) & + + (1.0_RKIND - (REAL(im)+1.0_RKIND)/M) * layerThicknessCur(:,iCell) ) & + config_fb_weight_3 & - * ( (im/M) * layerThicknessNew(:,iCell) & - + (1 - (im/M)) * layerThicknessCur(:,iCell) ) + * ( (REAL(im)/M) * layerThicknessNew(:,iCell) & + + (1.0_RKIND - REAL(im)/M) * layerThicknessCur(:,iCell) ) end do call mpas_timer_stop('FB_LTS interface prediction') From 683147cebedaa7bbae1ba59acc4160c2e1544852 Mon Sep 17 00:00:00 2001 From: Jeremy Lilly Date: Mon, 20 Nov 2023 09:08:13 -0800 Subject: [PATCH 24/34] Add implicit vmix step at the end of FB_LTS --- .../mpas_ocn_time_integration_fblts.F | 24 ++++++++++++++++++- 1 file changed, 23 insertions(+), 1 deletion(-) diff --git a/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_fblts.F b/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_fblts.F index 772caa6b19e7..8aed94b9002d 100644 --- a/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_fblts.F +++ b/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_fblts.F @@ -1321,7 +1321,29 @@ subroutine ocn_time_integrator_fblts(domain, dt)!{{{ call mpas_timer_stop("FB_LTS thickness halo update") call mpas_timer_stop("FB_LTS prognostic halo update") - + + ! Implicit vmix, DO WE NEED THIS? + call mpas_timer_start("FB_LTS implicit vmix") + if (.not. config_disable_vel_vmix) then + call ocn_diagnostic_solve(dt, statePool, forcingPool, meshPool, & + verticalMeshPool, scratchPool, tracersPool, 2) + call ocn_vmix_implicit(dt, meshPool, statePool, forcingPool, scratchPool, err, 2) + + ! Halo update + call mpas_timer_start("FB_LTS prognostic halo update") + + call mpas_timer_start("FB_LTS velocity halo update") + call mpas_dmpar_field_halo_exch(domain, 'normalVelocity', timeLevel=2) + call mpas_timer_stop("FB_LTS velocity halo update") + + call mpas_timer_start("FB_LTS thickness halo update") + call mpas_dmpar_field_halo_exch(domain, 'layerThickness', timeLevel=2) + call mpas_timer_stop("FB_LTS thickness halo update") + + call mpas_timer_stop("FB_LTS prognostic halo update") + end if + call mpas_timer_stop("FB_LTS implicit vmix") + call mpas_timer_stop("FB_LTS main loop") !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! From b2215c32691a94dbbfe67e100f6a54f168ca1ef6 Mon Sep 17 00:00:00 2001 From: Jeremy Lilly Date: Mon, 20 Nov 2023 14:05:58 -0800 Subject: [PATCH 25/34] Only advance interface 1 during the final step of coarse advancement --- .../mpas_ocn_time_integration_fblts.F | 26 ++++++++++--------- 1 file changed, 14 insertions(+), 12 deletions(-) diff --git a/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_fblts.F b/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_fblts.F index 8aed94b9002d..c55c09e26272 100644 --- a/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_fblts.F +++ b/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_fblts.F @@ -36,7 +36,7 @@ module ocn_time_integration_fblts use ocn_mesh use ocn_vmix use ocn_config - !use ocn_diagnostics_variables + use ocn_diagnostics_variables use ocn_equation_of_state use ocn_time_average_coupled use ocn_time_varying_forcing @@ -698,7 +698,7 @@ subroutine ocn_time_integrator_fblts(domain, dt)!{{{ do ic = 1, nCellsInLTSRegion(1,3) iCell = cellsInLTSRegion(1,3,ic) layerThicknessNew(:,iCell) = layerThicknessCur(:,iCell) & - + dt * layerThicknessTend(:,iCell) + + dt * layerThicknessTend(:,iCell) end do call mpas_timer_stop("FB_LTS advance thickness solution") @@ -750,13 +750,11 @@ subroutine ocn_time_integrator_fblts(domain, dt)!{{{ normalVelocityTend(:,:) = normalVelocityTend(:,:) & + normalVelocityTendSlow(:,:) - ! Interface regions - do iRegion = 1,nRegions - do ie = 1, nEdgesInLTSRegion(iRegion,2) - iEdge = edgesInLTSRegion(iRegion,2,ie) - normalVelocityNew(:,iEdge) = normalVelocityCur(:,iEdge) & - + dt * normalVelocityTend(:,iEdge) - end do + ! Interface 1 + do ie = 1, nEdgesInLTSRegion(1,2) + iEdge = edgesInLTSRegion(1,2,ie) + normalVelocityNew(:,iEdge) = normalVelocityCur(:,iEdge) & + + dt * normalVelocityTend(:,iEdge) end do ! Coarse do ie = 1, nEdgesInLTSRegion(2,1) @@ -1274,6 +1272,8 @@ subroutine ocn_time_integrator_fblts(domain, dt)!{{{ ! END: Normal velocity stage 3 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! COPY NEW INTO CUR ON FINE + end do !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! END FINE ADVANCEMENT @@ -1308,8 +1308,9 @@ subroutine ocn_time_integrator_fblts(domain, dt)!{{{ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! END INTERFACE CORRECTION !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + - ! Final prognostic halo update + ! Final prognostic halo update, is this needed? call mpas_timer_start("FB_LTS prognostic halo update") call mpas_timer_start("FB_LTS velocity halo update") @@ -1389,6 +1390,7 @@ subroutine ocn_time_integrator_fblts(domain, dt)!{{{ SSHGradient(indexSSHGradientMeridional, iCell) = gradSSHMeridional(iCell) end do + ! NEED TO GET CORRECT SSH AT TL==2!!???? call ocn_time_average_coupled_accumulate(statePool, forcingPool, 2) do iCell = 1, nCellsAll @@ -1624,7 +1626,7 @@ subroutine ocn_lts_thick_tend(LTSPool, tendPool, & computeOnInterfaceTwo, & computeOnCoarse - real (kind=RKIND), dimension(:,:), pointer :: & + real (kind=RKIND), dimension(:,:), pointer, intent(in) :: & layerThickness, & normalVelocity @@ -1795,7 +1797,7 @@ subroutine ocn_lts_vel_tend(LTSPool, tendPool, & computeOnInterfaceTwo, & computeOnCoarse - real (kind=RKIND), dimension(:,:), pointer :: & + real (kind=RKIND), dimension(:,:), pointer, intent(in) :: & layerThickness, & normalVelocity From 79b56897e48dfcc440c6b5949f1ed6b4ee0ed8f9 Mon Sep 17 00:00:00 2001 From: Jeremy Lilly Date: Mon, 20 Nov 2023 16:02:33 -0800 Subject: [PATCH 26/34] Copy new solution into current solution at the end of each fine substep, also calculate ssh before a call to diagnostics during the cleanup phase --- .../mpas_ocn_time_integration_fblts.F | 52 +++++++++++++++++-- 1 file changed, 49 insertions(+), 3 deletions(-) diff --git a/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_fblts.F b/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_fblts.F index c55c09e26272..37e81e3b7cd9 100644 --- a/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_fblts.F +++ b/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_fblts.F @@ -151,6 +151,10 @@ subroutine ocn_time_integrator_fblts(domain, dt)!{{{ layerThicknessSecondStage, & ! layer thickness at second stage of LTS layerThicknessForTend ! extra variable to store averages of data for tends calculations + ! Local pointer for ssh + real (kind=RKIND), dimension(:), pointer :: & + ssh + ! LTS objects real (kind=RKIND) :: & dtFine, & ! fine dt, defined as dt / M @@ -1272,7 +1276,37 @@ subroutine ocn_time_integrator_fblts(domain, dt)!{{{ ! END: Normal velocity stage 3 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - ! COPY NEW INTO CUR ON FINE + + ! Copy new solution into cur solution for use at the beginning + ! of the next substep + ! We might be able to speed this up by only doing a copy for + ! layer thickness and changing New to Cur in the above + ! stage 3 step for the normal velocity + call mpas_timer_start("FB_LTS copy new soln into cur soln") + + ! Fine adjacent to interface 1 + do ie = 1, nEdgesInLTSRegion(1,3) + iEdge = edgesInLTSRegion(1,3,ie) + normalVelocityCur(:,iEdge) = normalVelocityNew(:,iEdge) + end do + ! Fine in the interior, away from interface 1 + do ie = 1, nEdgesInLTSRegion(1,1) + iEdge = edgesInLTSRegion(1,1,ie) + normalVelocityCur(:,iEdge) = normalVelocityNew(:,iEdge) + end do + + ! Fine adjacent to interface 1 + do ic = 1, nCellsInLTSRegion(1,3) + iCell = cellsInLTSRegion(1,3,ic) + layerThicknessCur(:,iCell) = layerThicknessNew(:,iCell) + end do + ! Fine in the interior, away from interface 1 + do ic = 1, nCellsInLTSRegion(1,1) + iCell = cellsInLTSRegion(1,1,ic) + layerThicknessCur(:,iCell) = layerThicknessNew(:,iCell) + end do + + call mpas_timer_stop("FB_LTS copy new soln into cur soln") end do !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! @@ -1372,6 +1406,20 @@ subroutine ocn_time_integrator_fblts(domain, dt)!{{{ normalTransportVelocity(:,iEdge) = normalVelocityNew(:,iEdge) end do + ! Calculate ssh at new time level (2) before computing the diagnostics + ! It seems like we need this, as the call to ocn_diagnostic_solve + ! looks for ssh at time level 2 without calculating it itself + call mpas_pool_get_array(statePool, 'ssh', ssh, 2) + do iCell = 1, nCellsAll + k = maxLevelCell(iCell) + zTop(k:nVertLevels,iCell) = -bottomDepth(iCell) + layerThicknessNew(k,iCell) + do k = maxLevelCell(iCell)-1, minLevelCell(iCell), -1 + zTop(k,iCell) = zTop(k+1,iCell) + layerThicknessNew(k,iCell) + end do + ! copy zTop(1,iCell) into sea-surface height array + ssh(iCell) = zTop(minLevelCell(iCell),iCell) + end do + call mpas_reconstruct(meshPool, normalVelocityNew, & velocityX, velocityY, velocityZ, & velocityZonal, velocityMeridional, & @@ -1390,7 +1438,6 @@ subroutine ocn_time_integrator_fblts(domain, dt)!{{{ SSHGradient(indexSSHGradientMeridional, iCell) = gradSSHMeridional(iCell) end do - ! NEED TO GET CORRECT SSH AT TL==2!!???? call ocn_time_average_coupled_accumulate(statePool, forcingPool, 2) do iCell = 1, nCellsAll @@ -1401,7 +1448,6 @@ subroutine ocn_time_integrator_fblts(domain, dt)!{{{ SSHGradient(indexSSHGradientMeridional, iCell) = gradSSHMeridional(iCell) end do - ! Update diagnostics call ocn_diagnostic_solve(dt, statePool, forcingPool, meshPool, & verticalMeshPool, scratchPool, tracersPool, 2) From 85d25696d29a107824787c8d820c4aba305a53a0 Mon Sep 17 00:00:00 2001 From: Jeremy Lilly Date: Mon, 20 Nov 2023 16:27:06 -0800 Subject: [PATCH 27/34] Add FB_LTS to conidtionals in some tend calculations to match with LTS --- components/mpas-ocean/src/shared/mpas_ocn_diagnostics.F | 6 ++++-- components/mpas-ocean/src/shared/mpas_ocn_tendency.F | 4 +++- .../mpas-ocean/src/shared/mpas_ocn_vel_hadv_coriolis.F | 4 ++++ .../mpas-ocean/src/shared/mpas_ocn_vel_pressure_grad.F | 3 ++- .../mpas-ocean/src/shared/mpas_ocn_vel_tidal_potential.F | 3 ++- 5 files changed, 15 insertions(+), 5 deletions(-) diff --git a/components/mpas-ocean/src/shared/mpas_ocn_diagnostics.F b/components/mpas-ocean/src/shared/mpas_ocn_diagnostics.F index 5d3df864d25e..8f6905fc4e4f 100644 --- a/components/mpas-ocean/src/shared/mpas_ocn_diagnostics.F +++ b/components/mpas-ocean/src/shared/mpas_ocn_diagnostics.F @@ -4464,8 +4464,10 @@ subroutine ocn_diagnostics_init(domain, err)!{{{ ke_cell_flag = 1 endif - if ((trim(config_time_integrator) == 'RK4') .or.(trim(config_time_integrator) == 'LTS') ) then - ! For RK4 or LTS, PV includes f: PV = (eta+f)/h. + if ( (trim(config_time_integrator) == 'RK4') & + .or. (trim(config_time_integrator) == 'LTS') & + .or. (trim(config_time_integrator) == 'FB_LTS') ) then + ! For RK4, LTS, or FB_LTS, PV includes f: PV = (eta+f)/h. fCoef = 1 elseif (trim(config_time_integrator) == 'split_explicit' & .or.trim(config_time_integrator) == 'unsplit_explicit' & diff --git a/components/mpas-ocean/src/shared/mpas_ocn_tendency.F b/components/mpas-ocean/src/shared/mpas_ocn_tendency.F index a806c57412ef..860035d2600e 100644 --- a/components/mpas-ocean/src/shared/mpas_ocn_tendency.F +++ b/components/mpas-ocean/src/shared/mpas_ocn_tendency.F @@ -437,7 +437,9 @@ subroutine ocn_tend_vel(domain, tendPool, statePool, forcingPool, & ! Add tidal potential (if needed) call ocn_compute_tidal_potential_forcing(err) - if ((config_time_integrator == 'RK4') .or. (config_time_integrator =='LTS')) then + if ( (config_time_integrator == 'RK4') & + .or. (config_time_integrator =='LTS') & + .or. (config_time_integrator == 'FB_LTS') ) then ! for split explicit, tidal forcing is added in barotropic subcycles call ocn_vel_tidal_potential_tend(ssh,surfacePressure, tendVel, err) endif diff --git a/components/mpas-ocean/src/shared/mpas_ocn_vel_hadv_coriolis.F b/components/mpas-ocean/src/shared/mpas_ocn_vel_hadv_coriolis.F index a6cf77c9c527..0f4ad95757a3 100644 --- a/components/mpas-ocean/src/shared/mpas_ocn_vel_hadv_coriolis.F +++ b/components/mpas-ocean/src/shared/mpas_ocn_vel_hadv_coriolis.F @@ -359,6 +359,10 @@ subroutine ocn_vel_hadv_coriolis_init(err)!{{{ case ('LTS','lts') ! For LTS, coriolis tendency term includes f: (eta+f)/h. usePlanetVorticity = .true. + + case ('FB_LTS','fb_lts') + ! For FB_LTS, coriolis tendency term includes f: (eta+f)/h. + usePlanetVorticity = .true. case ('split_explicit') ! For split explicit, Coriolis tendency uses eta/h because diff --git a/components/mpas-ocean/src/shared/mpas_ocn_vel_pressure_grad.F b/components/mpas-ocean/src/shared/mpas_ocn_vel_pressure_grad.F index 0615856f6673..8e23913960f5 100644 --- a/components/mpas-ocean/src/shared/mpas_ocn_vel_pressure_grad.F +++ b/components/mpas-ocean/src/shared/mpas_ocn_vel_pressure_grad.F @@ -709,7 +709,8 @@ subroutine ocn_vel_pressure_grad_init(err)!{{{ end select - if (config_time_integrator == 'LTS') then + if ( (config_time_integrator == 'LTS') & + .or. (config_time_integrator == 'FB_LTS') ) then timeIntegratorLTS = .true. else timeIntegratorLTS = .false. diff --git a/components/mpas-ocean/src/shared/mpas_ocn_vel_tidal_potential.F b/components/mpas-ocean/src/shared/mpas_ocn_vel_tidal_potential.F index a13105200ed5..092c52c54d09 100644 --- a/components/mpas-ocean/src/shared/mpas_ocn_vel_tidal_potential.F +++ b/components/mpas-ocean/src/shared/mpas_ocn_vel_tidal_potential.F @@ -503,7 +503,8 @@ subroutine ocn_vel_tidal_potential_init(domain,err)!{{{ call mpas_log_write(' ') end do - if (config_time_integrator == 'LTS') then + if ( (config_time_integrator == 'LTS') & + .or. (config_time_integrator == 'FB_LTS') ) then timeIntegratorLTS = .true. else timeIntegratorLTS = .false. From 418fe141979113167b5cb821f61539f8d74407df Mon Sep 17 00:00:00 2001 From: Jeremy Lilly Date: Wed, 22 Nov 2023 10:08:12 -0800 Subject: [PATCH 28/34] Move addition of the slow tends to the fast tends to inside the tend calculation timer, also some general cleanup --- .../mpas_ocn_time_integration_fblts.F | 131 ++++-------------- 1 file changed, 30 insertions(+), 101 deletions(-) diff --git a/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_fblts.F b/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_fblts.F index 37e81e3b7cd9..7d48ceda5457 100644 --- a/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_fblts.F +++ b/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_fblts.F @@ -11,9 +11,9 @@ ! !> \brief MPAS barotropic ocean LTS Time integration scheme !> \author Jeremy Lilly -!> \date November 2022 +!> \date October 2023 !> \details -!> This module contains the LTS init routine and the LTS +!> This module contains the FB_LTS init routine and the FB_LTS !> barotropic ocean time integration scheme with splitting !> on the fast and slow tendency terms. ! @@ -72,7 +72,7 @@ module ocn_time_integration_fblts subroutine ocn_time_integrator_fblts(domain, dt)!{{{ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! Advance model state forward in time by the specified time step - ! using a local time stepping scheme with spltting of the fast and + ! using a local time stepping scheme with splitting of the fast and ! slow tendencies ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! @@ -118,12 +118,8 @@ subroutine ocn_time_integrator_fblts(domain, dt)!{{{ type (mpas_pool_type), pointer :: & LTSPool, & ! structure holding LTS variables tendSlowPool, & ! structure holding the slow tendency variables - tendSum1stPool, & ! structure holding one of the correction terms for the interface - tendSum2ndPool, & ! structure holding one of the correction terms for the interface tendSum3rdPool, & ! structure holding one of the correction terms for the interface prevTendSlowPool, nextTendSlowPool, & ! structures containing intermediate data - prevTendSum1stPool, nextTendSum1stPool, & ! structures containing intermediate data - prevTendSum2ndPool, nextTendSum2ndPool, & ! structures containing intermediate data prevTendSum3rdPool, nextTendSum3rdPool ! structures containing intermediate data ! Tend Array Pointers @@ -131,10 +127,6 @@ subroutine ocn_time_integrator_fblts(domain, dt)!{{{ normalVelocityTend, & ! normal velocity fast tendency layerThicknessTend, & ! layer thickness tendency normalVelocityTendSlow, & ! normal velocity slow tendency - normalVelocityTendSum1st, & ! one of the normal velocity correction terms for the interface - layerThicknessTendSum1st, & ! one of the layer thickness correction terms for the interface - normalVelocityTendSum2nd, & ! one of the normal velocity correction terms for the interface - layerThicknessTendSum2nd, & ! one of the layer thickness correction terms for the interface normalVelocityTendSum3rd, & ! one of the normal velocity correction terms for the interface layerThicknessTendSum3rd ! one of the layer thickness correction terms for the interface @@ -157,8 +149,7 @@ subroutine ocn_time_integrator_fblts(domain, dt)!{{{ ! LTS objects real (kind=RKIND) :: & - dtFine, & ! fine dt, defined as dt / M - alpha, alphaHat, beta, betaHat, gam, gamHat ! interpolation coefficients for the predictor step of LTS + dtFine ! fine dt, defined as dt / M integer, dimension(:,:), pointer :: & nCellsInLTSRegion, & ! number of cells in a given LTS region nEdgesInLTSRegion ! number of edges in a given LTS region @@ -166,9 +157,9 @@ subroutine ocn_time_integrator_fblts(domain, dt)!{{{ cellsInLTSRegion, & ! list of cells in a given LTS region edgesInLTSRegion ! list of edges in a given LTS region real (kind=RKIND) :: & - weight1st, weight2nd, weight3rd ! coefficents each RK stage + weight1st, weight2nd ! coefficients for each RK stage real (kind=RKIND) :: & - weightTendSum1st, weightTendSum2nd, weightTendSum3rd ! coefficients for the interface correction + weightTendSum3rd ! coefficients for the interface correction integer err err = 0 @@ -191,10 +182,9 @@ subroutine ocn_time_integrator_fblts(domain, dt)!{{{ nRegions = 2 dtFine = dt / M - ! Weights for RK stages + ! Weights for RK stages, weight for third stage is 1 weight1st = 1.0_RKIND / 3.0_RKIND weight2nd = 1.0_RKIND / 2.0_RKIND - weight3rd = 1.0_RKIND block => domain % blocklist @@ -241,39 +231,21 @@ subroutine ocn_time_integrator_fblts(domain, dt)!{{{ call mpas_pool_get_array(LTSPool, 'nEdgesInLTSRegion', nEdgesInLTSRegion) ! Create and retrieve additional pools for LTS - !call mpas_pool_create_pool(tendSum1stPool) - !call mpas_pool_clone_pool(tendPool, tendSum1stPool, 1) - !call mpas_pool_create_pool(tendSum2ndPool) - !call mpas_pool_clone_pool(tendPool, tendSum2ndPool, 1) call mpas_pool_create_pool(tendSum3rdPool) call mpas_pool_clone_pool(tendPool, tendSum3rdPool, 1) call mpas_pool_create_pool(tendSlowPool) call mpas_pool_clone_pool(tendPool, tendSlowPool, 1) - !call mpas_pool_add_subpool(block % structs, 'tend_sum_1st', tendSum1stPool) - !call mpas_pool_add_subpool(block % structs, 'tend_sum_2nd', tendSum2ndPool) call mpas_pool_add_subpool(block % structs, 'tend_sum_3rd', tendSum3rdPool) call mpas_pool_add_subpool(block % structs, 'tend_slow', tendSlowPool) call mpas_pool_get_array(tendSlowPool, 'normalVelocity', & normalVelocityTendSlow) - !call mpas_pool_get_array(tendSum1stPool, 'normalVelocity', & - ! normalVelocityTendSum1st) - !call mpas_pool_get_array(tendSum1stPool, 'layerThickness', & - ! layerThicknessTendSum1st) - !call mpas_pool_get_array(tendSum2ndPool, 'normalVelocity', & - ! normalVelocityTendSum2nd) - !call mpas_pool_get_array(tendSum2ndPool, 'layerThickness', & - ! layerThicknessTendSum2nd) call mpas_pool_get_array(tendSum3rdPool, 'normalVelocity', & normalVelocityTendSum3rd) call mpas_pool_get_array(tendSum3rdPool, 'layerThickness', & layerThicknessTendSum3rd) - ! not sure if we need this yet -- also need to dealloc if used! - !allocate(normalVelocityForTend(nVertLevels,nEdgesAll)) - !allocate(layerThicknessForTend(nVertLevels,nCellsAll)) - ! Init variables do iEdge = 1, nEdgesAll ! do we need this? do k = 1, maxLevelEdgeTop(iEdge) @@ -291,64 +263,28 @@ subroutine ocn_time_integrator_fblts(domain, dt)!{{{ end do end do - !normalVelocityTendSum1st(:,:) = 0.0_RKIND - !layerThicknessTendSum1st(:,:) = 0.0_RKIND - - !normalVelocityTendSum2nd(:,:) = 0.0_RKIND - !layerThicknessTendSum2nd(:,:) = 0.0_RKIND - normalVelocityTendSum3rd(:,:) = 0.0_RKIND layerThicknessTendSum3rd(:,:) = 0.0_RKIND if (associated(block % prev)) then - !call mpas_pool_get_subpool(block % prev % structs, 'tend_sum_1st', tendSum1stPool) - !call mpas_pool_get_subpool(block % prev % structs, 'tend_sum_2nd', tendSum2ndPool) call mpas_pool_get_subpool(block % prev % structs, 'tend_sum_3rd', tendSum3rdPool) call mpas_pool_get_subpool(block % prev % structs, 'tend_slow', tendSlowPool) else - !nullify(prevTendSum1stPool) - !nullify(prevTendSum2ndPool) nullify(prevTendSum3rdPool) nullify(prevTendSlowPool) end if if (associated(block % next)) then - !call mpas_pool_get_subpool(block % next % structs, 'tend_sum_1st', nextTendSum1stPool) - !call mpas_pool_get_subpool(block % next % structs, 'tend_sum_2nd', nextTendSum2ndPool) call mpas_pool_get_subpool(block % next % structs, 'tend_sum_3rd', nextTendSum3rdPool) call mpas_pool_get_subpool(block % next % structs, 'tend_slow', nextTendSlowPool) else - !nullify(nextTendSum1stPool) - !nullify(nextTendSum2ndPool) nullify(nextTendSum3rdPool) nullify(nextTendSlowPool) end if - !call mpas_pool_get_subpool(block % structs, 'tend_sum_1st', tendSum1stPool) - !call mpas_pool_get_subpool(block % structs, 'tend_sum_2nd', tendSum2ndPool) call mpas_pool_get_subpool(block % structs, 'tend_sum_3rd', tendSum3rdPool) call mpas_pool_get_subpool(block % structs, 'tend_slow', tendSlowPool) - !if (associated(prevTendSum1stPool) .and. associated(nextTendSum1stPool)) then - ! call mpas_pool_link_pools(tendSum1stPool, prevTendSum1stPool, nextTendSum1stPool) - !else if (associated(prevTendSum1stPool)) then - ! call mpas_pool_link_pools(tendSum1stPool, prevTendSum1stPool) - !else if (associated(nextTendSum1stPool)) then - ! call mpas_pool_link_pools(tendSum1stPool,nextPool=nextTendSum1stPool) - !else - ! call mpas_pool_link_pools(tendSum1stPool) - !end if - - !if (associated(prevTendSum2ndPool) .and. associated(nextTendSum2ndPool)) then - ! call mpas_pool_link_pools(tendSum2ndPool, prevTendSum2ndPool, nextTendSum2ndPool) - !else if (associated(prevTendSum2ndPool)) then - ! call mpas_pool_link_pools(tendSum2ndPool, prevTendSum2ndPool) - !else if (associated(nextTendSum2ndPool)) then - ! call mpas_pool_link_pools(tendSum2ndPool,nextPool=nextTendSum2ndPool) - !else - ! call mpas_pool_link_pools(tendSum2ndPool) - !end if - if (associated(prevTendSum3rdPool) .and. associated(nextTendSum3rdPool)) then call mpas_pool_link_pools(tendSum3rdPool, prevTendSum3rdPool, nextTendSum3rdPool) else if (associated(prevTendSum3rdPool)) then @@ -369,8 +305,6 @@ subroutine ocn_time_integrator_fblts(domain, dt)!{{{ call mpas_pool_link_pools(tendSlowPool) end if - !call mpas_pool_link_parinfo(block, tendSum1stPool) - !call mpas_pool_link_parinfo(block, tendSum2ndPool) call mpas_pool_link_parinfo(block, tendSum3rdPool) call mpas_pool_link_parinfo(block, tendSlowPool) @@ -483,6 +417,9 @@ subroutine ocn_time_integrator_fblts(domain, dt)!{{{ normalVelocityCur, layerThicknessForTend, & 0, 1, 1, 1, 1) + normalVelocityTend(:,:) = normalVelocityTend(:,:) & + + normalVelocityTendSlow(:,:) + call mpas_timer_stop("FB_LTS compute fast velocity tendencies") call mpas_timer_stop("FB_LTS compute fast tendencies") @@ -490,9 +427,6 @@ subroutine ocn_time_integrator_fblts(domain, dt)!{{{ call mpas_timer_start("FB_LTS advance solution") call mpas_timer_start("FB_LTS advance velocity solution") - normalVelocityTend(:,:) = normalVelocityTend(:,:) & - + normalVelocityTendSlow(:,:) - ! Interface regions do iRegion = 1,nRegions do ie = 1, nEdgesInLTSRegion(iRegion,2) @@ -613,15 +547,15 @@ subroutine ocn_time_integrator_fblts(domain, dt)!{{{ normalVelocityFirstStage, layerThicknessForTend, & 0, 1, 1, 1, 1) + normalVelocityTend(:,:) = normalVelocityTend(:,:) & + + normalVelocityTendSlow(:,:) + call mpas_timer_stop("FB_LTS compute fast velocity tendencies") call mpas_timer_stop("FB_LTS compute fast tendencies") ! Advance normal velocity solution with first stage call mpas_timer_start("FB_LTS advance solution") call mpas_timer_start("FB_LTS advance velocity solution") - - normalVelocityTend(:,:) = normalVelocityTend(:,:) & - + normalVelocityTendSlow(:,:) ! Interface regions do iRegion = 1,nRegions @@ -744,15 +678,15 @@ subroutine ocn_time_integrator_fblts(domain, dt)!{{{ normalVelocitySecondStage, layerThicknessForTend, & 0, 0, 1, 0, 1) + normalVelocityTend(:,:) = normalVelocityTend(:,:) & + + normalVelocityTendSlow(:,:) + call mpas_timer_stop("FB_LTS compute fast velocity tendencies") call mpas_timer_stop("FB_LTS compute fast tendencies") ! Advance normal velocity solution with first stage call mpas_timer_start("FB_LTS advance solution") call mpas_timer_start("FB_LTS advance velocity solution") - - normalVelocityTend(:,:) = normalVelocityTend(:,:) & - + normalVelocityTendSlow(:,:) ! Interface 1 do ie = 1, nEdgesInLTSRegion(1,2) @@ -916,15 +850,15 @@ subroutine ocn_time_integrator_fblts(domain, dt)!{{{ normalVelocityForTend, layerThicknessForTend, & 1, 1, 0, 0, 0) + normalVelocityTend(:,:) = normalVelocityTend(:,:) & + + normalVelocityTendSlow(:,:) + call mpas_timer_stop("FB_LTS compute fast velocity tendencies") call mpas_timer_stop("FB_LTS compute fast tendencies") ! Advance normal velocity solution with first stage call mpas_timer_start("FB_LTS advance solution") call mpas_timer_start("FB_LTS advance velocity solution") - - normalVelocityTend(:,:) = normalVelocityTend(:,:) & - + normalVelocityTendSlow(:,:) ! Fine cells adjacent to interface 1 do ie = 1, nEdgesInLTSRegion(1,3) @@ -1070,6 +1004,9 @@ subroutine ocn_time_integrator_fblts(domain, dt)!{{{ normalVelocityForTend, layerThicknessForTend, & 1, 1, 0, 0, 0) + normalVelocityTend(:,:) = normalVelocityTend(:,:) & + + normalVelocityTendSlow(:,:) + call mpas_timer_stop("FB_LTS compute fast velocity tendencies") call mpas_timer_stop("FB_LTS compute fast tendencies") @@ -1077,9 +1014,6 @@ subroutine ocn_time_integrator_fblts(domain, dt)!{{{ call mpas_timer_start("FB_LTS advance solution") call mpas_timer_start("FB_LTS advance velocity solution") - normalVelocityTend(:,:) = normalVelocityTend(:,:) & - + normalVelocityTendSlow(:,:) - ! Fine cells adjacent to interface 1 do ie = 1, nEdgesInLTSRegion(1,3) iEdge = edgesInLTSRegion(1,3,ie) @@ -1239,6 +1173,9 @@ subroutine ocn_time_integrator_fblts(domain, dt)!{{{ normalVelocityForTend, layerThicknessForTend, & 1, 1, 1, 1, 0) + normalVelocityTend(:,:) = normalVelocityTend(:,:) & + + normalVelocityTendSlow(:,:) + call mpas_timer_stop("FB_LTS compute fast velocity tendencies") call mpas_timer_stop("FB_LTS compute fast tendencies") @@ -1246,9 +1183,6 @@ subroutine ocn_time_integrator_fblts(domain, dt)!{{{ call mpas_timer_start("FB_LTS advance solution") call mpas_timer_start("FB_LTS advance velocity solution") - normalVelocityTend(:,:) = normalVelocityTend(:,:) & - + normalVelocityTendSlow(:,:) - ! Increment interface correction terms do iRegion = 1,nRegions do ie = 1, nEdgesInLTSRegion(iRegion,2) @@ -1452,13 +1386,9 @@ subroutine ocn_time_integrator_fblts(domain, dt)!{{{ call ocn_diagnostic_solve(dt, statePool, forcingPool, meshPool, & verticalMeshPool, scratchPool, tracersPool, 2) - !call mpas_pool_destroy_pool(tendSum1stPool) - !call mpas_pool_destroy_pool(tendSum2ndPool) call mpas_pool_destroy_pool(tendSum3rdPool) call mpas_pool_destroy_pool(tendSlowPool) - !call mpas_pool_remove_subpool(block % structs, 'tend_sum_1st') - !call mpas_pool_remove_subpool(block % structs, 'tend_sum_2nd') call mpas_pool_remove_subpool(block % structs, 'tend_sum_3rd') call mpas_pool_remove_subpool(block % structs, 'tend_slow') @@ -1473,7 +1403,7 @@ end subroutine ocn_time_integrator_fblts!}}} ! !> \brief Initialization for MPAS ocean FB_LTS time integration scheme !> \author Giacomo Capodaglio -!> \date October 2023 +!> \date November 2022 !> \details !> This routine computes the LTS arrays ! @@ -1651,7 +1581,8 @@ end subroutine ocn_time_integration_fblts_init!}}} !> \author Jeremy Lilly !> \date October 2023 !> \details -!> This routine calculates the tendencies on different LTS regions +!> This routine calculates the thickness tendency on different +!> LTS regions ! !----------------------------------------------------------------------- @@ -1822,7 +1753,8 @@ end subroutine ocn_lts_thick_tend!}}} !> \author Jeremy Lilly !> \date October 2023 !> \details -!> This routine calculates the tendencies on different LTS regions +!> This routine calculates the velocity tendency on different +!> LTS regions ! !----------------------------------------------------------------------- @@ -1891,9 +1823,6 @@ subroutine ocn_lts_vel_tend(LTSPool, tendPool, & call mpas_pool_get_array(tendPool, 'normalVelocity', normalVelocityTend) - !call mpas_pool_get_array(statePool, 'ssh', & - ! ssh, timeLevelVelocity) - if (config_use_self_attraction_loading) then ssh_sal_on = 1.0_RKIND else From 4c6ef8da5150a150bebb1b73705a20bd151b0c0c Mon Sep 17 00:00:00 2001 From: Jeremy Lilly Date: Wed, 22 Nov 2023 10:41:35 -0800 Subject: [PATCH 29/34] Remove uneeded modules --- .../src/mode_forward/mpas_ocn_time_integration_fblts.F | 7 ------- 1 file changed, 7 deletions(-) diff --git a/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_fblts.F b/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_fblts.F index 7d48ceda5457..8de7a482be01 100644 --- a/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_fblts.F +++ b/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_fblts.F @@ -21,25 +21,18 @@ module ocn_time_integration_fblts - use mpas_derived_types use mpas_pool_routines - use mpas_constants use mpas_dmpar use mpas_threading use mpas_vector_reconstruction - use mpas_spline_interpolation use mpas_timer - use ocn_constants use ocn_tendency use ocn_diagnostics use ocn_mesh use ocn_vmix use ocn_config - use ocn_diagnostics_variables - use ocn_equation_of_state use ocn_time_average_coupled - use ocn_time_varying_forcing implicit none private From d153f427b3527521022aa9cd28672d4e35b7defe Mon Sep 17 00:00:00 2001 From: Jeremy Lilly Date: Wed, 22 Nov 2023 11:12:25 -0800 Subject: [PATCH 30/34] Remove unneeded halo update --- .../src/mode_forward/mpas_ocn_time_integration_fblts.F | 10 ---------- 1 file changed, 10 deletions(-) diff --git a/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_fblts.F b/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_fblts.F index 8de7a482be01..df7f6b0be891 100644 --- a/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_fblts.F +++ b/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_fblts.F @@ -697,16 +697,6 @@ subroutine ocn_time_integrator_fblts(domain, dt)!{{{ call mpas_timer_stop("FB_LTS advance velocity solution") call mpas_timer_stop("FB_LTS advance solution") - ! Halo update - ! Might be able to remove this? - call mpas_timer_start("FB_LTS prognostic halo update") - - call mpas_timer_start("FB_LTS velocity halo update") - call mpas_dmpar_field_halo_exch(domain, 'normalVelocity', timeLevel=2) - call mpas_timer_stop("FB_LTS velocity halo update") - - call mpas_timer_stop("FB_LTS prognostic halo update") - ! END: Normal velocity stage 3 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! From cf01384a52d1ecaed556286aa59dd9cd8e423d0c Mon Sep 17 00:00:00 2001 From: Jeremy Lilly Date: Mon, 4 Dec 2023 10:25:09 -0800 Subject: [PATCH 31/34] Remove unneeded calculation of ssh during cleanup phase --- .../mode_forward/mpas_ocn_time_integration_fblts.F | 14 -------------- 1 file changed, 14 deletions(-) diff --git a/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_fblts.F b/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_fblts.F index df7f6b0be891..ee483014ea3a 100644 --- a/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_fblts.F +++ b/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_fblts.F @@ -1323,20 +1323,6 @@ subroutine ocn_time_integrator_fblts(domain, dt)!{{{ normalTransportVelocity(:,iEdge) = normalVelocityNew(:,iEdge) end do - ! Calculate ssh at new time level (2) before computing the diagnostics - ! It seems like we need this, as the call to ocn_diagnostic_solve - ! looks for ssh at time level 2 without calculating it itself - call mpas_pool_get_array(statePool, 'ssh', ssh, 2) - do iCell = 1, nCellsAll - k = maxLevelCell(iCell) - zTop(k:nVertLevels,iCell) = -bottomDepth(iCell) + layerThicknessNew(k,iCell) - do k = maxLevelCell(iCell)-1, minLevelCell(iCell), -1 - zTop(k,iCell) = zTop(k+1,iCell) + layerThicknessNew(k,iCell) - end do - ! copy zTop(1,iCell) into sea-surface height array - ssh(iCell) = zTop(minLevelCell(iCell),iCell) - end do - call mpas_reconstruct(meshPool, normalVelocityNew, & velocityX, velocityY, velocityZ, & velocityZonal, velocityMeridional, & From b617c9c30bcbe26ace3fd17d54c361ed1756d2a1 Mon Sep 17 00:00:00 2001 From: Jeremy Lilly Date: Mon, 4 Dec 2023 10:45:44 -0800 Subject: [PATCH 32/34] Rearrange how the copying of the solutions between New and Cur is done at the end of each subcycle, for performance --- .../mpas_ocn_time_integration_fblts.F | 48 ++++++++++++------- 1 file changed, 30 insertions(+), 18 deletions(-) diff --git a/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_fblts.F b/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_fblts.F index ee483014ea3a..33a4c0bc3285 100644 --- a/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_fblts.F +++ b/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_fblts.F @@ -1163,6 +1163,10 @@ subroutine ocn_time_integrator_fblts(domain, dt)!{{{ call mpas_timer_stop("FB_LTS compute fast tendencies") ! Advance normal velocity solution with third stage + ! Here, we write into normalVelocityCur, rather than + ! normalVelocityNew, so that the New data is stored + ! in normalVelocityCur to be used at the beginning of the + ! next subcycled time-step call mpas_timer_start("FB_LTS advance solution") call mpas_timer_start("FB_LTS advance velocity solution") @@ -1177,13 +1181,13 @@ subroutine ocn_time_integrator_fblts(domain, dt)!{{{ ! Fine cells adjacent to interface 1 do ie = 1, nEdgesInLTSRegion(1,3) iEdge = edgesInLTSRegion(1,3,ie) - normalVelocityNew(:,iEdge) = normalVelocityCur(:,iEdge) & + normalVelocityCur(:,iEdge) = normalVelocityCur(:,iEdge) & + dtFine * normalVelocityTend(:,iEdge) end do ! Fine in the interior, away from interface 1 do ie = 1, nEdgesInLTSRegion(1,1) iEdge = edgesInLTSRegion(1,1,ie) - normalVelocityNew(:,iEdge) = normalVelocityCur(:,iEdge) & + normalVelocityCur(:,iEdge) = normalVelocityCur(:,iEdge) & + dtFine * normalVelocityTend(:,iEdge) end do @@ -1196,22 +1200,11 @@ subroutine ocn_time_integrator_fblts(domain, dt)!{{{ ! Copy new solution into cur solution for use at the beginning ! of the next substep - ! We might be able to speed this up by only doing a copy for - ! layer thickness and changing New to Cur in the above - ! stage 3 step for the normal velocity - call mpas_timer_start("FB_LTS copy new soln into cur soln") + ! We only do this for layerThickness here because we write + ! into normalVelocityCur rather than normalVelocityNew above + ! for stage 3 of the subcycled normal velocity + call mpas_timer_start("FB_LTS copy soln") - ! Fine adjacent to interface 1 - do ie = 1, nEdgesInLTSRegion(1,3) - iEdge = edgesInLTSRegion(1,3,ie) - normalVelocityCur(:,iEdge) = normalVelocityNew(:,iEdge) - end do - ! Fine in the interior, away from interface 1 - do ie = 1, nEdgesInLTSRegion(1,1) - iEdge = edgesInLTSRegion(1,1,ie) - normalVelocityCur(:,iEdge) = normalVelocityNew(:,iEdge) - end do - ! Fine adjacent to interface 1 do ic = 1, nCellsInLTSRegion(1,3) iCell = cellsInLTSRegion(1,3,ic) @@ -1223,9 +1216,28 @@ subroutine ocn_time_integrator_fblts(domain, dt)!{{{ layerThicknessCur(:,iCell) = layerThicknessNew(:,iCell) end do - call mpas_timer_stop("FB_LTS copy new soln into cur soln") + call mpas_timer_stop("FB_LTS copy soln") + end do ! end of fine time-step subcycling loop + + + ! Copy data from normalVelocityCur into normalVelocityNew + ! which currently holds the New data. + call mpas_timer_start("FB_LTS copy soln") + + ! Fine adjacent to interface 1 + do ie = 1, nEdgesInLTSRegion(1,3) + iEdge = edgesInLTSRegion(1,3,ie) + normalVelocityNew(:,iEdge) = normalVelocityCur(:,iEdge) + end do + ! Fine in the interior, away from interface 1 + do ie = 1, nEdgesInLTSRegion(1,1) + iEdge = edgesInLTSRegion(1,1,ie) + normalVelocityNew(:,iEdge) = normalVelocityCur(:,iEdge) end do + + call mpas_timer_stop("FB_LTS copy soln") + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! END FINE ADVANCEMENT !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! From b2de77325fa42f630aec09cd6c011ebf1f255f5e Mon Sep 17 00:00:00 2001 From: Jeremy Lilly Date: Mon, 4 Dec 2023 14:07:15 -0800 Subject: [PATCH 33/34] Add upwind option for thickness flux --- .../mpas_ocn_time_integration_fblts.F | 65 ++++++++++++++----- 1 file changed, 49 insertions(+), 16 deletions(-) diff --git a/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_fblts.F b/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_fblts.F index 33a4c0bc3285..6db62099a9f3 100644 --- a/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_fblts.F +++ b/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_fblts.F @@ -1627,24 +1627,57 @@ subroutine ocn_lts_thick_tend(LTSPool, tendPool, & call mpas_pool_get_array(tendPool, 'layerThickness', layerThicknessTend) ! calculate layerThickEdgeFlux - ! add upwind formulation, ask Giacomo about branch - do iEdge = 1, nEdgesAll - kmin = minLevelEdgeBot(iEdge) - kmax = maxLevelEdgeTop(iEdge) - cell1 = cellsOnEdge(1,iEdge) - cell2 = cellsOnEdge(2,iEdge) - do k = 1,nVertLevels - ! initialize layerThicknessEdgeMean to avoid divide by - ! zero and NaN problems. - layerThickEdgeFlux(k,iEdge) = -1.0e34_RKIND + if (config_thickness_flux_type == 'centered') then + + do iEdge = 1, nEdgesAll + kmin = minLevelEdgeBot(iEdge) + kmax = maxLevelEdgeTop(iEdge) + cell1 = cellsOnEdge(1,iEdge) + cell2 = cellsOnEdge(2,iEdge) + do k = 1,nVertLevels + ! initialize layerThicknessEdgeMean to avoid divide by + ! zero and NaN problems. + layerThickEdgeFlux(k,iEdge) = -1.0e34_RKIND + end do + do k = kmin,kmax + ! central differenced + layerThickEdgeFlux(k,iEdge) = 0.5_RKIND * & + ( layerThickness(k,cell1) + & + layerThickness(k,cell2) ) + end do end do - do k = kmin,kmax - ! central differenced - layerThickEdgeFlux(k,iEdge) = 0.5_RKIND * & - ( layerThickness(k,cell1) & - + layerThickness(k,cell2) ) + + else if (config_thickness_flux_type == 'upwind') then + + do iEdge = 1, nEdgesAll + kmin = minLevelEdgeBot(iEdge) + kmax = maxLevelEdgeTop(iEdge) + cell1 = cellsOnEdge(1,iEdge) + cell2 = cellsOnEdge(2,iEdge) + do k=1,nVertLevels + ! initialize layerThicknessEdgeFlux to avoid divide by + ! zero and NaN problems. + layerThickEdgeFlux(k,iEdge) = -1.0e34_RKIND + end do + do k = kmin,kmax + if (normalVelocity(k,iEdge) > 0.0_RKIND) then + layerThickEdgeFlux(k,iEdge) = layerThickness(k,cell1) + elseif (normalVelocity(k,iEdge) < 0.0_RKIND) then + layerThickEdgeFlux(k,iEdge) = layerThickness(k,cell2) + else + layerThickEdgeFlux(k,iEdge) = max(layerThickness(k,cell1), & + layerThickness(k,cell2)) + end if + end do end do - end do + + else + + call mpas_log_write('ERROR: config_thickness_flux_type selected not implemented for LTS', & + messageType=MPAS_LOG_CRIT) + call abort + + end if layerThicknessTend(:, :) = 0.0_RKIND From afcac9aad624d43df6bb2fb2d81e3879655da645 Mon Sep 17 00:00:00 2001 From: Jeremy Lilly Date: Mon, 4 Dec 2023 14:36:05 -0800 Subject: [PATCH 34/34] Check whether tidal potential forcing is on before adding the beta SAL term to the ssh gradient --- .../mpas_ocn_time_integration_fblts.F | 29 +++++++++++-------- 1 file changed, 17 insertions(+), 12 deletions(-) diff --git a/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_fblts.F b/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_fblts.F index 6db62099a9f3..36b6dc05766f 100644 --- a/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_fblts.F +++ b/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_fblts.F @@ -1673,7 +1673,7 @@ subroutine ocn_lts_thick_tend(LTSPool, tendPool, & else - call mpas_log_write('ERROR: config_thickness_flux_type selected not implemented for LTS', & + call mpas_log_write('ERROR: config_thickness_flux_type selected not implemented for FB_LTS', & messageType=MPAS_LOG_CRIT) call abort @@ -1824,7 +1824,7 @@ subroutine ocn_lts_vel_tend(LTSPool, tendPool, & ic, i, iCell, kmin, kmax real (kind=RKIND) :: & - invdcEdge, betaSelfAttrLoad, flux, ssh_sal_on + invdcEdge, betaSelfAttrLoad, flux, ssh_sal_on, tidal_pot_for_on !----------------------------------------------------------------- ! Begin routine @@ -1842,6 +1842,11 @@ subroutine ocn_lts_vel_tend(LTSPool, tendPool, & else ssh_sal_on = 0.0_RKIND endif + if (config_use_tidal_potential_forcing) then + tidal_pot_for_on = 1.0_RKIND + else + tidal_pot_for_on = 0.0_RKIND + end if ! ssh do iCell = 1, nCellsAll @@ -1869,8 +1874,8 @@ subroutine ocn_lts_vel_tend(LTSPool, tendPool, & normalVelocityTend(k,iEdge) = normalVelocityTend(k,iEdge) & - edgeMask(k,iEdge) * invdcEdge & * ( gravity * ( (ssh(cell2) - ssh(cell1)) & - - (1.0_RKIND - ssh_sal_on) * betaSelfAttrLoad & - * (ssh(cell2) - ssh(cell1)) ) ) + - tidal_pot_for_on * (1.0_RKIND - ssh_sal_on) & + * betaSelfAttrLoad * (ssh(cell2) - ssh(cell1)) ) ) end do end do end if @@ -1888,8 +1893,8 @@ subroutine ocn_lts_vel_tend(LTSPool, tendPool, & normalVelocityTend(k,iEdge) = normalVelocityTend(k,iEdge) & - edgeMask(k,iEdge) * invdcEdge & * ( gravity * ( (ssh(cell2) - ssh(cell1)) & - - (1.0_RKIND - ssh_sal_on) * betaSelfAttrLoad & - * (ssh(cell2) - ssh(cell1)) ) ) + - tidal_pot_for_on * (1.0_RKIND - ssh_sal_on) & + * betaSelfAttrLoad * (ssh(cell2) - ssh(cell1)) ) ) end do end do end if @@ -1907,8 +1912,8 @@ subroutine ocn_lts_vel_tend(LTSPool, tendPool, & normalVelocityTend(k,iEdge) = normalVelocityTend(k,iEdge) & - edgeMask(k,iEdge) * invdcEdge & * ( gravity * ( (ssh(cell2) - ssh(cell1)) & - - (1.0_RKIND - ssh_sal_on) * betaSelfAttrLoad & - * (ssh(cell2) - ssh(cell1)) ) ) + - tidal_pot_for_on * (1.0_RKIND - ssh_sal_on) & + * betaSelfAttrLoad * (ssh(cell2) - ssh(cell1)) ) ) end do end do end if @@ -1926,8 +1931,8 @@ subroutine ocn_lts_vel_tend(LTSPool, tendPool, & normalVelocityTend(k,iEdge) = normalVelocityTend(k,iEdge) & - edgeMask(k,iEdge) * invdcEdge & * ( gravity * ( (ssh(cell2) - ssh(cell1)) & - - (1.0_RKIND - ssh_sal_on) * betaSelfAttrLoad & - * (ssh(cell2) - ssh(cell1)) ) ) + - tidal_pot_for_on * (1.0_RKIND - ssh_sal_on) & + * betaSelfAttrLoad * (ssh(cell2) - ssh(cell1)) ) ) end do end do end if @@ -1945,8 +1950,8 @@ subroutine ocn_lts_vel_tend(LTSPool, tendPool, & normalVelocityTend(k,iEdge) = normalVelocityTend(k,iEdge) & - edgeMask(k,iEdge) * invdcEdge & * ( gravity * ( (ssh(cell2) - ssh(cell1)) & - - (1.0_RKIND - ssh_sal_on) * betaSelfAttrLoad & - * (ssh(cell2) - ssh(cell1)) ) ) + - tidal_pot_for_on * (1.0_RKIND - ssh_sal_on) & + * betaSelfAttrLoad * (ssh(cell2) - ssh(cell1)) ) ) end do end do end if