diff --git a/components/mpas-ocean/bld/build-namelist b/components/mpas-ocean/bld/build-namelist index 5fee738461f8..f749d4df0043 100755 --- a/components/mpas-ocean/bld/build-namelist +++ b/components/mpas-ocean/bld/build-namelist @@ -511,6 +511,7 @@ add_default($nl, 'config_transport_tests_flow_id'); add_default($nl, 'config_dt'); add_default($nl, 'config_time_integrator'); +add_default($nl, 'config_triskCV'); ######################## # Namelist group: hmix # diff --git a/components/mpas-ocean/bld/build-namelist-section b/components/mpas-ocean/bld/build-namelist-section index d7d4efe740f0..0352598fd83f 100644 --- a/components/mpas-ocean/bld/build-namelist-section +++ b/components/mpas-ocean/bld/build-namelist-section @@ -54,6 +54,7 @@ add_default($nl, 'config_transport_tests_flow_id'); add_default($nl, 'config_dt'); add_default($nl, 'config_time_integrator'); +add_default($nl, 'config_triskCV'); ######################## # Namelist group: hmix # diff --git a/components/mpas-ocean/bld/namelist_files/namelist_defaults_mpaso.xml b/components/mpas-ocean/bld/namelist_files/namelist_defaults_mpaso.xml index 789acecc6e54..f45b1b70e027 100644 --- a/components/mpas-ocean/bld/namelist_files/namelist_defaults_mpaso.xml +++ b/components/mpas-ocean/bld/namelist_files/namelist_defaults_mpaso.xml @@ -53,6 +53,7 @@ '00:10:00' '00:30:00' 'split_explicit' +.true. .false. diff --git a/components/mpas-ocean/bld/namelist_files/namelist_definition_mpaso.xml b/components/mpas-ocean/bld/namelist_files/namelist_definition_mpaso.xml index 6e0153ead725..573cad3461c8 100644 --- a/components/mpas-ocean/bld/namelist_files/namelist_definition_mpaso.xml +++ b/components/mpas-ocean/bld/namelist_files/namelist_definition_mpaso.xml @@ -235,6 +235,13 @@ Valid values: 'split_explicit', 'RK4', 'unsplit_explicit', 'split_implicit' Default: Defined in namelist_defaults.xml + +Turns on the computation of the kinetic energy in the TRiSK-CV way. + +Valid values: .true. or .false. +Default: Defined in namelist_defaults.xml + diff --git a/components/mpas-ocean/src/Registry.xml b/components/mpas-ocean/src/Registry.xml index 9e9128f39dda..08a16297cbf4 100644 --- a/components/mpas-ocean/src/Registry.xml +++ b/components/mpas-ocean/src/Registry.xml @@ -208,6 +208,10 @@ description="Time integration method." possible_values="'split_explicit', 'RK4', 'unsplit_explicit', 'split_implicit'" /> + + + + + + + + - + + + + + @@ -2260,6 +2290,9 @@ + @@ -2304,8 +2337,8 @@ /> - + + + 1e-10) then + if (config_apvm_scale_factor>1e-10 .or. config_aust_scale_factor>1e-10) then allocate(vorticityGradientNormalComponent(nVertLevels, nEdges), & vorticityGradientTangentialComponent(nVertLevels, nEdges)) @@ -1106,15 +1111,17 @@ subroutine ocn_diagnostic_solve_vorticity(dt, normalVelocity, tangentialVelocity !$acc loop seq #endif do k = minLevelEdgeTop(iEdge), maxLevelEdgeBot(iEdge) - vorticityGradientTangentialComponent(k,iEdge) = & - (normalizedRelativeVorticityVertex(k,vertex2) - normalizedRelativeVorticityVertex(k,vertex1)) * invLength + vorticityGradientTangentialComponent(k,iEdge) = & + (normalizedRelativeVorticityVertex(k,vertex2) - normalizedRelativeVorticityVertex(k,vertex1)) * invLength enddo enddo #ifndef MPAS_OPENACC !$omp end do + !$omp end parallel #endif + if (config_apvm_scale_factor>1e-10) then ! ! Modify PV edge with upstream bias. ! @@ -1125,21 +1132,53 @@ subroutine ocn_diagnostic_solve_vorticity(dt, normalVelocity, tangentialVelocity #else !$omp do schedule(runtime) private(k) #endif - do iEdge = 1,nEdges + do iEdge = 1,nEdges +#ifdef MPAS_OPENACC + !$acc loop seq +#endif + do k = minLevelEdgeTop(iEdge), maxLevelEdgeBot(iEdge) + normalizedRelativeVorticityEdge(k,iEdge) = normalizedRelativeVorticityEdge(k,iEdge) & + - apvm_scale_factor * dt * & + ( normalVelocity(k,iEdge) * vorticityGradientNormalComponent(k,iEdge) & + + tangentialVelocity(k,iEdge) * vorticityGradientTangentialComponent(k,iEdge) ) + end do + enddo +#ifndef MPAS_OPENACC + !$omp end do + !$omp end parallel +#endif + end if + + if (config_aust_scale_factor>1e-10) then + ! + ! Modify PV edge with scaling given by cell width. + ! +#ifdef MPAS_OPENACC + !$acc parallel loop, & + !$acc present(minLevelEdgeTop, maxLevelEdgeBot, tangentialVelocity, & + !$acc normalVelocity, normalizedRelativeVorticityEdge) +#else + !$omp do schedule(runtime) private(k) +#endif + do iEdge = 1,nEdges #ifdef MPAS_OPENACC !$acc loop seq #endif - do k = minLevelEdgeTop(iEdge), maxLevelEdgeBot(iEdge) - normalizedRelativeVorticityEdge(k,iEdge) = normalizedRelativeVorticityEdge(k,iEdge) & - - apvm_scale_factor * dt * & - ( normalVelocity(k,iEdge) * vorticityGradientNormalComponent(k,iEdge) & - + tangentialVelocity(k,iEdge) * vorticityGradientTangentialComponent(k,iEdge) ) - end do - enddo + se = 0.5_RKIND * sqrt(dvEdge(iEdge)*dcEdge(iEdge)) + do k = minLevelEdgeTop(iEdge), maxLevelEdgeBot(iEdge) + invNormVel = 1.0_RKIND/sqrt(normalVelocity(k,iEdge)*normalVelocity(k,iEdge) + & + tangentialVelocity(k,iEdge)*tangentialVelocity(k,iEdge) + 1e-12) + normalizedRelativeVorticityEdge(k,iEdge) = normalizedRelativeVorticityEdge(k,iEdge) & + - aust_scale_factor * se * invNormVel * & + ( normalVelocity(k,iEdge) * vorticityGradientNormalComponent(k,iEdge) & + + tangentialVelocity(k,iEdge) * vorticityGradientTangentialComponent(k,iEdge) ) + end do + enddo #ifndef MPAS_OPENACC !$omp end do !$omp end parallel #endif + end if deallocate(vorticityGradientNormalComponent, & vorticityGradientTangentialComponent) @@ -1479,7 +1518,7 @@ subroutine ocn_diagnostic_solve_vortVel(relativeVorticity, & layerThicknessEdgeFlux, normalVelocity, normalTransportVelocity, & normalGMBolusVelocity, vertVelocityTop, vertTransportVelocityTop, & vertGMBolusVelocityTop, relativeVorticityCell, divergence, & - kineticEnergyCell, tangentialVelocity)!{{{ + kineticEnergyCell, kineticEnergyEdge, u_perp, tangentialVelocity)!{{{ implicit none @@ -1518,6 +1557,10 @@ subroutine ocn_diagnostic_solve_vortVel(relativeVorticity, & divergence real (kind=RKIND), dimension(:,:), intent(out) :: & kineticEnergyCell + real (kind=RKIND), dimension(:,:), intent(out) :: & + kineticEnergyEdge + real (kind=RKIND), dimension(:,:), intent(out) :: & + u_perp real (kind=RKIND), dimension(:,:), intent(out) :: & tangentialVelocity @@ -1528,13 +1571,17 @@ subroutine ocn_diagnostic_solve_vortVel(relativeVorticity, & !----------------------------------------------------------------- integer :: nCells, nEdges - integer :: iCell, iVertex, iEdge + integer :: iCell, iVertex, iEdge, jEdge integer :: i, j, k integer :: edgeSignOnCell_temp, eoe real(kind=RKIND) :: invAreaCell1 - real(kind=RKIND) :: dcEdge_temp, dvEdge_temp, r_tmp, weightsOnEdge_temp + real(kind=RKIND) :: dcEdge_temp, dvEdge_temp, r_tmp, weightsOnEdge_temp, & + weightsOnEdge_lsqr_temp real (kind=RKIND), dimension(:), allocatable:: div_hu,div_huTransport,div_huGMBolus + + !call mpas_pool_get_array(diagnosticsPool, 'kineticEnergyCell', kineticEnergyCell, timeLevel) + #ifdef MPAS_OPENACC !$acc declare device_resident(div_hu, div_huTransport, div_huGMBolus) #endif @@ -1573,6 +1620,7 @@ subroutine ocn_diagnostic_solve_vortVel(relativeVorticity, & ! Need 0,1,2 halo cells nCells = nCellsHalo( 2 ) + nEdges = nEdgesHalo( 2 ) !nEdgesAll ! ! Compute divergence, kinetic energy, and vertical velocity @@ -1597,7 +1645,7 @@ subroutine ocn_diagnostic_solve_vortVel(relativeVorticity, & #endif do iCell = 1, nCells divergence(:, iCell) = 0.0_RKIND - kineticEnergyCell(:, iCell) = 0.0_RKIND + !kineticEnergyCell(:, iCell) = 0.0_RKIND div_hu(:) = 0.0_RKIND div_huTransport(:) = 0.0_RKIND div_huGMBolus(:) = 0.0_RKIND @@ -1612,8 +1660,8 @@ subroutine ocn_diagnostic_solve_vortVel(relativeVorticity, & divergence(k, iCell) = divergence(k, iCell) - edgeSignOnCell_temp * r_tmp div_hu(k) = div_hu(k) - layerThicknessEdgeFlux(k, iEdge) * edgeSignOnCell_temp * r_tmp - kineticEnergyCell(k, iCell) = kineticEnergyCell(k, iCell) & - + 0.25 * r_tmp * dcEdge_temp * normalVelocity(k,iEdge) + !kineticEnergyCell(k, iCell) = kineticEnergyCell(k, iCell) & + ! + 0.25 * r_tmp * dcEdge_temp * normalVelocity(k,iEdge) ! Compute vertical velocity from the horizontal total transport div_huTransport(k) = div_huTransport(k) & - layerThicknessEdgeFlux(k, iEdge) * edgeSignOnCell_temp & @@ -1637,35 +1685,173 @@ subroutine ocn_diagnostic_solve_vortVel(relativeVorticity, & !$omp end do !$omp end parallel #endif - - deallocate(div_hu,div_huTransport,div_huGMBolus) - - nEdges = nEdgesHalo( 2 ) +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + if (config_triskCV) then #ifdef MPAS_OPENACC !$acc parallel loop gang vector & - !$acc present(tangentialVelocity, nEdgesOnEdge, edgesOnEdge, weightsOnEdge, & + !$acc present(tangentialVelocity, nEdgesOnEdge, edgesOnEdge, & + !$acc weightsOnEdge, weightsOnEdge_lsqr, & !$acc maxLevelEdgeTop, normalVelocity, minLevelEdgeBot) & !$acc private(i, eoe, weightsOnEdge_temp, k) #else !$omp parallel !$omp do schedule(runtime) private(i, eoe, weightsOnEdge_temp, k) #endif - do iEdge = 1, nEdges - tangentialVelocity(:, iEdge) = 0.0_RKIND - ! Compute v (tangential) velocities - do i = 1, nEdgesOnEdge(iEdge) - eoe = edgesOnEdge(i,iEdge) - weightsOnEdge_temp = weightsOnEdge(i, iEdge) - do k = minLevelEdgeBot(iEdge), maxLevelEdgeTop(iEdge) - tangentialVelocity(k,iEdge) = tangentialVelocity(k,iEdge) & + do iEdge = 1, nEdges + tangentialVelocity(:, iEdge) = 0.0_RKIND + ! Compute v (tangential) velocities + do i = 1, nEdgesOnEdge(iEdge) + eoe = edgesOnEdge(i,iEdge) + weightsOnEdge_temp = weightsOnEdge_lsqr(i, iEdge) + do k = minLevelEdgeBot(iEdge), maxLevelEdgeTop(iEdge) + tangentialVelocity(k,iEdge) = tangentialVelocity(k,iEdge) & + weightsOnEdge_temp * normalVelocity(k, eoe) + end do + end do + end do +#ifndef MPAS_OPENACC + !$omp end do + !$omp end parallel +#endif + else +#ifdef MPAS_OPENACC + !$acc parallel loop gang vector & + !$acc present(tangentialVelocity, nEdgesOnEdge, edgesOnEdge, & + !$acc weightsOnEdge, weightsOnEdge_lsqr, & + !$acc maxLevelEdgeTop, normalVelocity, minLevelEdgeBot) & + !$acc private(i, eoe, weightsOnEdge_temp, k) +#else + !$omp parallel + !$omp do schedule(runtime) private(i, eoe, weightsOnEdge_temp, k) +#endif + do iEdge = 1, nEdges + tangentialVelocity(:, iEdge) = 0.0_RKIND + ! Compute v (tangential) velocities + do i = 1, nEdgesOnEdge(iEdge) + eoe = edgesOnEdge(i,iEdge) + weightsOnEdge_temp = weightsOnEdge(i, iEdge) + do k = minLevelEdgeBot(iEdge), maxLevelEdgeTop(iEdge) + tangentialVelocity(k,iEdge) = tangentialVelocity(k,iEdge) & + + weightsOnEdge_temp * normalVelocity(k, eoe) + end do end do end do - end do #ifndef MPAS_OPENACC !$omp end do !$omp end parallel #endif + end if + + if (config_triskCV) then + kineticEnergyEdge(:,:) = 0.0_RKIND + +!#ifdef MPAS_OPENACC +! !$acc parallel loop gang vector & +! !$acc present(u_perp, nEdgesOnEdge, edgesOnEdge, weightsOnEdge_lsqr, & +! !$acc minLevelEdgeBot, maxLevelEdgeTop, normalVelocity) +! !$acc private(iEdge, nEdges, j, jEdge, weightsOnEdge_lsqr_temp, k) +!#else +! !$omp parallel +! !$omp do schedule(runtime) private(iEdge, nEdges, j, jEdge, weightsOnEdge_lsqr_temp, k) +!#endif +! do iEdge = 1, nEdges +! u_perp(:,iEdge) = 0.0_RKIND +! do j = 1, nEdgesOnEdge(iEdge) +! jEdge = edgesOnEdge(j,iEdge) +! weightsOnEdge_lsqr_temp = weightsOnEdge_lsqr(j,iEdge) +! do k = minLevelEdgeBot(iEdge), maxLevelEdgeTop(iEdge) +! u_perp(k,iEdge) = u_perp(k,iEdge) + weightsOnEdge_lsqr_temp*normalVelocity(k,jEdge) +! end do +! end do +! end do +!#ifndef MPAS_OPENACC +! !$omp end do +! !$omp end parallel +!#endif + +#ifdef MPAS_OPENACC + !$acc parallel loop gang vector & + !$acc present(minLevelEdgeBot, maxLevelEdgeTop, kineticEnergyEdge, & + !$acc normalVelocity, u_perp) + !$acc private(iEdge, nEdges, k) +#else + !$omp parallel + !$omp do schedule(runtime) private(iEdge, nEdges, k) +#endif + do iEdge = 1, nEdges + do k = minLevelEdgeBot(iEdge), maxLevelEdgeTop(iEdge) + !kineticEnergyEdge(k,iEdge) = 0.5_RKIND*(normalVelocity(k,iEdge)*normalVelocity(k,iEdge) + u_perp(k,iEdge)*u_perp(k,iEdge)) + kineticEnergyEdge(k,iEdge) = 0.5_RKIND*(normalVelocity(k,iEdge)*normalVelocity(k,iEdge) + & + tangentialVelocity(k,iEdge)*tangentialVelocity(k,iEdge)) + end do + end do +#ifndef MPAS_OPENACC + !$omp end do + !$omp end parallel +#endif + +#ifdef MPAS_OPENACC + !$acc parallel loop gang vector & + !$acc present(kineticEnergyCell, nEdgesOnCell, edgesOnCell, & + !$acc minLevelCell, maxLevelCell, dcEdge, dvEdge, areaCell) + !$acc private(iCell, nCells, i, iEdge, k) +#else + !$omp parallel + !$omp do schedule(runtime) private(iCell, nCells, i, iEdge, k) +#endif + do iCell = 1, nCells + kineticEnergyCell(:, iCell) = 0.0_RKIND + do i = 1, nEdgesOnCell(iCell) + iEdge = edgesOnCell(i, iCell) + do k = minLevelCell(iCell), maxLevelCell(iCell) + kineticEnergyCell(k,iCell) = kineticEnergyCell(k,iCell) + 0.25_RKIND * dcEdge(iEdge) * dvEdge(iEdge) * kineticEnergyEdge(k,iEdge) + end do + end do + do k = minLevelCell(iCell), maxLevelCell(iCell) + kineticEnergyCell(k,iCell) = kineticEnergyCell(k,iCell) / areaCell(iCell) + end do + end do +#ifndef MPAS_OPENACC + !$omp end do + !$omp end parallel +#endif + + else + +#ifdef MPAS_OPENACC + !$acc parallel loop gang vector & + !$acc present(kineticEnergyCell, invAreaCell, nEdgesOnCell, & + !$acc edgesOnCell, dcEdge, dvEdge, minLevelCell, & + !$acc maxLevelCell, normalVelocity) + !$acc private(iCell, nCells, invAreaCell1, i, iEdge, dcEdge_temp, & + !$acc dvEdge_temp, r_tmp) +#else + !$omp parallel + !$omp do schedule(runtime) private(iCell, nCells, invAreaCell1, i, iEdge, & + !$omp dcEdge_temp, dvEdge_temp, r_tmp) +#endif + do iCell = 1, nCells + kineticEnergyCell(:, iCell) = 0.0_RKIND + invAreaCell1 = invAreaCell(iCell) + do i = 1, nEdgesOnCell(iCell) + iEdge = edgesOnCell(i, iCell) + dcEdge_temp = dcEdge(iEdge) + dvEdge_temp = dvEdge(iEdge) + do k = minLevelCell(iCell), maxLevelCell(iCell) + r_tmp = dvEdge_temp * normalVelocity(k, iEdge) * invAreaCell1 + kineticEnergyCell(k, iCell) = kineticEnergyCell(k, iCell) & + + 0.25 * r_tmp * dcEdge_temp * normalVelocity(k,iEdge) + end do + end do + end do +#ifndef MPAS_OPENACC + !$omp end do + !$omp end parallel +#endif + end if +!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + deallocate(div_hu,div_huTransport,div_huGMBolus) end subroutine ocn_diagnostic_solve_vortVel!}}} @@ -2375,6 +2561,19 @@ subroutine ocn_fuperp(statePool, meshPool, timeLevelIn)!{{{ normalVelocity(k,iEdge) = normalVelocity(k,iEdge) + weightsOnEdge(j,iEdge) * normalBaroclinicVelocity(k,eoe) & * fEdge(eoe) end do + + !if(k==5 .and. indexToEdgeID(iEdge)==1000) then + !print*, '---inside---' + ! do j = 1, nEdgesOnEdge(iEdge) + ! print*, indexToEdgeID(iEdge), k, j, normalVelocity(k,edgesOnEdge(j,iEdge)) + ! print*, "----associated normalBaroclinicVelocity" + ! do eoe = 1, nEdgesOnEdge(edgesOnEdge(j,iEdge)) + ! print*, eoe, normalBaroclinicVelocity(k,edgesOnEdge(eoe,edgesOnEdge(j,iEdge))) + ! end do + ! end do + !print*, '---end inside---' + !end if + end do end do #ifdef MPAS_OPENACC diff --git a/components/mpas-ocean/src/shared/mpas_ocn_diagnostics_variables.F b/components/mpas-ocean/src/shared/mpas_ocn_diagnostics_variables.F index 3325960bda3e..c8d5cb984bf0 100644 --- a/components/mpas-ocean/src/shared/mpas_ocn_diagnostics_variables.F +++ b/components/mpas-ocean/src/shared/mpas_ocn_diagnostics_variables.F @@ -59,6 +59,8 @@ module ocn_diagnostics_variables real (kind=RKIND), dimension(:,:), pointer :: layerThickEdgeFlux real (kind=RKIND), dimension(:,:), pointer :: layerThickEdgeMean real (kind=RKIND), dimension(:,:), pointer :: kineticEnergyCell + real (kind=RKIND), dimension(:,:), pointer :: kineticEnergyEdge + real (kind=RKIND), dimension(:,:), pointer :: u_perp real (kind=RKIND), dimension(:,:), pointer :: vertVelocityTop real (kind=RKIND), dimension(:,:), pointer :: vertTransportVelocityTop real (kind=RKIND), dimension(:,:), pointer :: vertGMBolusVelocityTop @@ -460,6 +462,10 @@ subroutine ocn_diagnostics_variables_init(domain, jenkinsOn, hollandJenkinsOn, e barotropicForcing) call mpas_pool_get_array(diagnosticsPool, 'kineticEnergyCell', & kineticEnergyCell) + call mpas_pool_get_array(diagnosticsPool, 'kineticEnergyEdge', & + kineticEnergyEdge) + call mpas_pool_get_array(diagnosticsPool, 'u_perp', & + u_perp) call mpas_pool_get_array(diagnosticsPool, 'transportVelocityX', & transportVelocityX) @@ -700,6 +706,8 @@ subroutine ocn_diagnostics_variables_init(domain, jenkinsOn, hollandJenkinsOn, e !$acc GMStreamFuncZonal, & !$acc relativeVorticity, & !$acc kineticEnergyCell, & + !$acc kineticEnergyEdge, & + !$acc u_perp, & !$acc normPlanetVortEdge, & !$acc velocityMeridional, & !$acc transportVelocityX, & @@ -940,6 +948,8 @@ subroutine ocn_diagnostics_variables_destroy(err) !{{{ !$acc GMStreamFuncZonal, & !$acc relativeVorticity, & !$acc kineticEnergyCell, & + !$acc kineticEnergyEdge, & + !$acc u_perp, & !$acc normPlanetVortEdge, & !$acc velocityMeridional, & !$acc transportVelocityX, & @@ -1138,6 +1148,8 @@ subroutine ocn_diagnostics_variables_destroy(err) !{{{ GMStreamFuncZonal, & relativeVorticity, & kineticEnergyCell, & + kineticEnergyEdge, & + u_perp, & normPlanetVortEdge, & velocityMeridional, & transportVelocityX, & diff --git a/components/mpas-ocean/src/shared/mpas_ocn_init_routines.F b/components/mpas-ocean/src/shared/mpas_ocn_init_routines.F index 65b839b1005b..9586767297c7 100644 --- a/components/mpas-ocean/src/shared/mpas_ocn_init_routines.F +++ b/components/mpas-ocean/src/shared/mpas_ocn_init_routines.F @@ -38,6 +38,7 @@ module ocn_init_routines use ocn_constants use ocn_config use ocn_mesh + use ocn_ke_triskcv use ocn_surface_land_ice_fluxes use ocn_forcing @@ -96,6 +97,7 @@ subroutine ocn_init_routines_block(block, dt, err)!{{{ real (kind=RKIND), intent(in) :: dt integer, intent(out) :: err + type (domain_type) :: domain type (mpas_pool_type), pointer :: meshPool type (mpas_pool_type), pointer :: statePool, tracersPool type (mpas_pool_type), pointer :: forcingPool, scratchPool @@ -169,7 +171,9 @@ subroutine ocn_init_routines_block(block, dt, err)!{{{ !$acc enter data copyin(landIceDraft) endif #endif + call ocn_diagnostic_solve(dt, statePool, forcingPool, meshPool, scratchPool, tracersPool) + #ifdef MPAS_OPENACC !$acc update host(layerThickEdgeFlux, layerThickEdgeMean) !$acc update host(relativeVorticity, circulation) @@ -209,6 +213,18 @@ subroutine ocn_init_routines_block(block, dt, err)!{{{ endif !$acc exit data delete(layerThickness, normalVelocity) #endif + + !if (config_triskCV) then + ! ! Compute the areas of the diamond-shaped control volumes needed in + ! ! TRiSK-CV + ! call controlVolumeEdge(domain, meshPool) + ! ! Compute a matrix needed in TRiSK-CV for the computation of the + ! ! kinetic energy + ! call buildInvMatrixBtB(domain, meshPool) + ! ! Compute the weightsOnEdge needed in TRiSK-CV for the computation of + ! ! the kinetic energy in TRiSK-CV + ! call Wlsqr(domain, meshPool) + !end if layerThickness(:, nCells+1) = 0.0_RKIND diff --git a/components/mpas-ocean/src/shared/mpas_ocn_ke_triskCV.F b/components/mpas-ocean/src/shared/mpas_ocn_ke_triskCV.F new file mode 100644 index 000000000000..2f2f602359dc --- /dev/null +++ b/components/mpas-ocean/src/shared/mpas_ocn_ke_triskCV.F @@ -0,0 +1,342 @@ +! Copyright (c) 2013, 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_ke_triskCV +! +!> \brief MPAS ocean driver for the kinetic energy in TRiSK-CV +!> \author Sara Calandrini +!> \date 17 September 2021 +!> \details +!> This module contains the routines for computing the kinetic +!> energy in TRiSK-CV. +! +!----------------------------------------------------------------------- + +module ocn_ke_triskCV + + use mpas_timer + use mpas_derived_types + use mpas_pool_routines + use mpas_constants + use mpas_threading + use mpas_vector_reconstruction + use mpas_stream_manager + use mpas_io_units + use mpas_dmpar + + use ocn_constants + use ocn_config + use ocn_thick_ale + use ocn_diagnostics_variables + use ocn_mesh + + implicit none + private + save + + !-------------------------------------------------------------------- + ! + ! Public parameters + ! + !-------------------------------------------------------------------- + + !-------------------------------------------------------------------- + ! + ! Public member functions + ! + !-------------------------------------------------------------------- + + public :: Wlsqr, & + inv_3x3, & + sphericalTriangleArea, & + sphere_distance + + !-------------------------------------------------------------------- + ! + ! Private module variables + ! + !-------------------------------------------------------------------- + +!*********************************************************************** + +contains + + subroutine inv_3x3(amat,adim,ainv,vdim, adet) + + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! inv3x3 computes the inverse of a 3x3 matrix + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + implicit none + + !------------------------------------------- arguments ! + integer, intent( in) :: adim + real (kind=RKIND), intent( in) :: amat(adim,*) + integer, intent( in) :: vdim + real (kind=RKIND), intent(out) :: ainv(vdim,*) + real (kind=RKIND), intent(out) :: adet + + !------------------------------------------- variables ! + real (kind=RKIND) :: & + aa2233,aa2332,aa2133,aa2331,aa2132, & + aa2231,aa1233,aa1332,aa1223,aa1322, & + aa1133,aa1331,aa1123,aa1321,aa1132, & + aa1231,aa1122,aa1221 + + !------------------------------------------- form A^-1 ! + aa2233 = amat(2,2) * amat(3,3) + aa2332 = amat(2,3) * amat(3,2) + aa2133 = amat(2,1) * amat(3,3) + aa2331 = amat(2,3) * amat(3,1) + aa2132 = amat(2,1) * amat(3,2) + aa2231 = amat(2,2) * amat(3,1) + + adet = & + amat(1,1) * (aa2233 - aa2332) - & + amat(1,2) * (aa2133 - aa2331) + & + amat(1,3) * (aa2132 - aa2231) + + aa1233 = amat(1,2) * amat(3,3) + aa1332 = amat(1,3) * amat(3,2) + aa1223 = amat(1,2) * amat(2,3) + aa1322 = amat(1,3) * amat(2,2) + aa1133 = amat(1,1) * amat(3,3) + aa1331 = amat(1,3) * amat(3,1) + aa1123 = amat(1,1) * amat(2,3) + aa1321 = amat(1,3) * amat(2,1) + aa1132 = amat(1,1) * amat(3,2) + aa1231 = amat(1,2) * amat(3,1) + aa1122 = amat(1,1) * amat(2,2) + aa1221 = amat(1,2) * amat(2,1) + + ainv(1,1) = (aa2233 - aa2332) + ainv(1,2) = -(aa1233 - aa1332) + ainv(1,3) = (aa1223 - aa1322) + + ainv(2,1) = -(aa2133 - aa2331) + ainv(2,2) = (aa1133 - aa1331) + ainv(2,3) = -(aa1123 - aa1321) + + ainv(3,1) = (aa2132 - aa2231) + ainv(3,2) = -(aa1132 - aa1231) + ainv(3,3) = (aa1122 - aa1221) + + return + + end subroutine inv_3x3 + + subroutine Wlsqr(domain) + + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! buildMatrixBtB computes the least square weightsOnEdge needed for the + ! computation of the kinetic energy in the TRiSK++ discretization + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + implicit none + + type (domain_type), intent(inout) :: domain + integer :: iEdge, jEdge, iEoE, i, j, k, cell1 + integer, dimension(2) :: iEdgeCells, jEdgeCells + real (kind=RKIND) :: re_x, re_y, re_z, ne_x, ne_y, ne_z, norm, xTangent, yTangent, zTangent, normTangent, detinv + real (kind=RKIND), dimension(:), pointer :: areaCVEdge + real (kind=RKIND), dimension(:,:,:), pointer :: invBtB + + real (kind=RKIND) :: tmat(25,3), matt(3,25), rmat(3,3), rinv(3,3), & + rdet, xmul, ymul, zmul + real (kind=RKIND), dimension(:), allocatable :: elen + real (kind=RKIND), dimension(:,:), allocatable :: enrm, edir, eprp + + allocate(enrm(nEdgesHalo(3),3)) + allocate(edir(nEdgesAll,3)) + allocate(eprp(nEdgesHalo(3),3)) + allocate(elen(nEdgesAll)) + +#ifdef MPAS_OPENACC + !$acc parallel loop gang vector & + !$acc present(nEdgesHalo, xVertex, yVertex, zVertex, verticesOnEdge, & + !$acc xEdge, yEdge, zEdge, sphereRadius, nEdgesOnEdge, & + !$acc edgesOnEdge, cellsOnEdge, nCellsAll, xCell, yCell, & + !$acc zCell, weightsOnEdge_lsqr, onSphere) & + !$acc private(iEdge, tmat, matt, rmat, rinv, eprp, jEdge, norm, & + !$acc enrm, iEoE, cell1, edir, elen, xmul, ymul, zmul, & + !$acc rdet) +#else + !$omp parallel + !$omp do schedule(runtime) & + !$omp private(iEdge, tmat, matt, rmat, rinv, eprp, jEdge, norm, & + !$omp enrm, iEoE, cell1, edir, elen, xmul, ymul, zmul, & + !$omp rdet) +#endif + !- build lsqr weights for each edge in the mesh + do iEdge = 1, nEdgesHalo(3) + tmat(:,:) = 0.0_RKIND + matt(:,:) = 0.0_RKIND + rmat(:,:) = 0.0_RKIND + rinv(:,:) = 0.0_RKIND + + !- unit tangent (on spherical surface) to each edge (in E^3) + eprp(iEdge,1) = xVertex(verticesOnEdge(2,iEdge)) - & + xVertex(verticesOnEdge(1,iEdge)) + eprp(iEdge,2) = yVertex(verticesOnEdge(2,iEdge)) - & + yVertex(verticesOnEdge(1,iEdge)) + eprp(iEdge,3) = zVertex(verticesOnEdge(2,iEdge)) - & + zVertex(verticesOnEdge(1,iEdge)) + + norm = sqrt(eprp(iEdge,1)**2 + eprp(iEdge,2)**2 + eprp(iEdge,3)**2) + + eprp(iEdge,1) = eprp(iEdge,1) / norm + eprp(iEdge,2) = eprp(iEdge,2) / norm + eprp(iEdge,3) = eprp(iEdge,3) / norm + + !- unit normals (to spherical surface) at each edge + if (onSphere) then + enrm(iEdge,1) = xEdge(iEdge) / sphereRadius + enrm(iEdge,2) = yEdge(iEdge) / sphereRadius + enrm(iEdge,3) = zEdge(iEdge) / sphereRadius + else + enrm(iEdge,1) = 0.0_RKIND + enrm(iEdge,2) = 0.0_RKIND + enrm(iEdge,3) = 1.0_RKIND + end if + + !- assemble lsqr matrix T: n_hat dot u = u_e + tmat(1,1) = enrm(iEdge,1) + tmat(1,2) = enrm(iEdge,2) + tmat(1,3) = enrm(iEdge,3) + + do jEdge = 1, nEdgesOnEdge(iEdge) + + iEoE = edgesOnEdge(jEdge,iEdge) + if(iEoE>nEdgesAll) cycle + if (cellsOnEdge(1,iEoE)>nCellsAll) then + edir(iEoE,1) = xCell(cellsOnEdge(2,iEoE)) - xEdge(iEoE) + edir(iEoE,2) = yCell(cellsOnEdge(2,iEoE)) - yEdge(iEoE) + edir(iEoE,3) = zCell(cellsOnEdge(2,iEoE)) - zEdge(iEoE) + elseif (cellsOnEdge(2,iEoE)>nCellsAll) then + edir(iEoE,1) = xEdge(iEoE) - xCell(cellsOnEdge(1,iEoE)) + edir(iEoE,2) = yEdge(iEoE) - yCell(cellsOnEdge(1,iEoE)) + edir(iEoE,3) = zEdge(iEoE) - zCell(cellsOnEdge(1,iEoE)) + else + edir(iEoE,1) = xCell(cellsOnEdge(2,iEoE)) - & + xCell(cellsOnEdge(1,iEoE)) + edir(iEoE,2) = yCell(cellsOnEdge(2,iEoE)) - & + yCell(cellsOnEdge(1,iEoE)) + edir(iEoE,3) = zCell(cellsOnEdge(2,iEoE)) - & + zCell(cellsOnEdge(1,iEoE)) + end if + + elen(iEoE) = sqrt(edir(iEoE,1)**2 + edir(iEoE,2)**2 + edir(iEoE,3)**2) + + edir(iEoE,1) = edir(iEoE,1) / elen(iEoE) + edir(iEoE,2) = edir(iEoE,2) / elen(iEoE) + edir(iEoE,3) = edir(iEoE,3) / elen(iEoE) + + tmat(jEdge+1,1) = edir(iEoE,1) + tmat(jEdge+1,2) = edir(iEoE,2) + tmat(jEdge+1,3) = edir(iEoE,3) + end do + + !- form lsqr matrices: (T'*T) * u = T' * u_e + matt = transpose(tmat) + + rmat = matmul(matt,tmat) + + !- factorise matrices: r_inv = (T'*T) ^ -1 + call inv_3x3(rmat, 3, rinv, 3, rdet) + + !- assemble lsqr weight: t_hat dot r_inv * T' + do jEdge = 1, nEdgesOnEdge(iEdge) + xmul = ( & + rinv(1,1) * matt(1,jEdge+1) + & + rinv(1,2) * matt(2,jEdge+1) + & + rinv(1,3) * matt(3,jEdge+1) & + ) / rdet + + ymul = ( & + rinv(2,1) * matt(1,jEdge+1) + & + rinv(2,2) * matt(2,jEdge+1) + & + rinv(2,3) * matt(3,jEdge+1) & + ) / rdet + + zmul = ( & + rinv(3,1) * matt(1,jEdge+1) + & + rinv(3,2) * matt(2,jEdge+1) + & + rinv(3,3) * matt(3,jEdge+1) & + ) / rdet + + weightsOnEdge_lsqr(jEdge,iEdge) = & + eprp(iEdge,1) * xmul + & + eprp(iEdge,2) * ymul + & + eprp(iEdge,3) * zmul + + end do + end do +#ifndef MPAS_OPENACC + !$omp end do + !$omp end parallel +#endif + + call mpas_dmpar_field_halo_exch(domain, 'weightsOnEdge_lsqr') + + deallocate(eprp) + deallocate(edir) + deallocate(enrm) + deallocate(elen) + + end subroutine Wlsqr + + + real function sphericalTriangleArea(laEdge, loEdge, laVertex, loVertex, laCell, loCell) + + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! sphericalTriangleArea uses the spherical analog of Heron's formula to + ! compute the area of a triangle on the surface of a sphere + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + implicit none + + real (kind=RKIND), intent(in) :: laEdge, loEdge, laVertex, loVertex, laCell, loCell + real (kind=RKIND) :: tanqe, s, a, b, c, one, radius + one = 1.0 + radius = 6371220.0 + + a = sphere_distance(laCell, loCell, laVertex, loVertex, one) + b = sphere_distance(laCell, loCell, laEdge, loEdge, one) + c = sphere_distance(laEdge, loEdge, laVertex, loVertex, one) + s = 0.5*(a+b+c) + + tanqe = sqrt(tan(0.5*s)*tan(0.5*(s-a))*tan(0.5*(s-b))*tan(0.5*(s-c))) + sphericalTriangleArea = 4.*atan(tanqe)*radius**2 + + end function sphericalTriangleArea + + real function sphere_distance(lat1, lon1, lat2, lon2, radius) + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! Compute the great-circle distance between (lat1, lon1) and (lat2, lon2) on + ! a sphere with given radius. + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + implicit none + + real (kind=RKIND), intent(in) :: lat1, lon1, lat2, lon2, radius + + real (kind=RKIND) :: arg1 + + arg1 = sqrt( sin(0.5*(lat2-lat1))**2 + & + cos(lat1)*cos(lat2)*sin(0.5*(lon2-lon1))**2 ) + sphere_distance = 2.*radius*asin(arg1) + + end function sphere_distance + + +end module ocn_ke_triskCV +!||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| +! vim: foldmethod=marker diff --git a/components/mpas-ocean/src/shared/mpas_ocn_mesh.F b/components/mpas-ocean/src/shared/mpas_ocn_mesh.F index 8b0a0e543fc3..e9fe904b8894 100644 --- a/components/mpas-ocean/src/shared/mpas_ocn_mesh.F +++ b/components/mpas-ocean/src/shared/mpas_ocn_mesh.F @@ -30,7 +30,7 @@ module ocn_mesh use mpas_derived_types use mpas_pool_routines use mpas_constants - use mpas_log + use mpas_log use ocn_config @@ -145,6 +145,7 @@ module ocn_mesh real(kind=RKIND), public, dimension(:, :), pointer :: & weightsOnEdge, &! weights on each edge + weightsOnEdge_lsqr, & kiteAreasOnVertex, &! real (vertexDegree nVertices) edgeTangentVectors, &! tangent unit vector at edge edgeNormalVectors, &! normal unit vector at edge @@ -437,6 +438,8 @@ subroutine ocn_meshCreate(domain) !{{{ areaTriangle) call mpas_pool_get_array(meshPool, 'weightsOnEdge', & weightsOnEdge) + call mpas_pool_get_array(meshPool, 'weightsOnEdge_lsqr', & + weightsOnEdge_lsqr) call mpas_pool_get_array(meshPool, 'bottomDepth', & bottomDepth) call mpas_pool_get_array(meshPool, 'refBottomDepth', & @@ -458,6 +461,8 @@ subroutine ocn_meshCreate(domain) !{{{ call mpas_pool_get_array(meshPool, 'weightsOnEdge', & weightsOnEdge) + call mpas_pool_get_array(meshPool, 'weightsOnEdge_lsqr', & + weightsOnEdge_lsqr) call mpas_pool_get_array(meshPool, 'kiteAreasOnVertex', & kiteAreasOnVertex) call mpas_pool_get_array(meshPool, 'edgeTangentVectors', & @@ -653,6 +658,7 @@ subroutine ocn_meshCreate(domain) !{{{ !$acc edgeSignOnCell, & !$acc edgeSignOnVertex, & !$acc weightsOnEdge, & + !$acc weightsOnEdge_lsqr, & !$acc kiteAreasOnVertex, & !$acc edgeTangentVectors, & !$acc edgeNormalVectors, & @@ -783,6 +789,7 @@ subroutine ocn_meshDestroy(err) !{{{ !$acc edgeSignOnCell, & !$acc edgeSignOnVertex, & !$acc weightsOnEdge, & + !$acc weightsOnEdge_lsqr, & !$acc kiteAreasOnVertex, & !$acc edgeTangentVectors,& !$acc edgeNormalVectors, & @@ -865,6 +872,7 @@ subroutine ocn_meshDestroy(err) !{{{ angleEdge, & distanceToCoast, & weightsOnEdge, & + weightsOnEdge_lsqr, & kiteAreasOnVertex, & edgeTangentVectors, & edgeNormalVectors, & @@ -1037,6 +1045,8 @@ subroutine ocn_meshUpdateFields(domain) !{{{ ! areaTriangle) !call mpas_pool_get_array(meshPool, 'weightsOnEdge', & ! weightsOnEdge) + !call mpas_pool_get_array(meshPool, 'weightsOnEdge_lsqr', & + ! weightsOnEdge_lsqr) !call mpas_pool_get_array(meshPool, 'bottomDepth', & ! bottomDepth) !call mpas_pool_get_array(meshPool, 'refBottomDepth', & @@ -1058,6 +1068,8 @@ subroutine ocn_meshUpdateFields(domain) !{{{ ! !call mpas_pool_get_array(meshPool, 'weightsOnEdge', & ! weightsOnEdge) + !call mpas_pool_get_array(meshPool, 'weightsOnEdge_lsqr', & + ! weightsOnEdge_lsqr) !call mpas_pool_get_array(meshPool, 'kiteAreasOnVertex', & ! kiteAreasOnVertex) !call mpas_pool_get_array(meshPool, 'edgeTangentVectors', & @@ -1214,6 +1226,7 @@ subroutine ocn_meshUpdateFields(domain) !{{{ !$acc edgeSignOnCell, & !$acc edgeSignOnVertex, & !$acc weightsOnEdge, & + !$acc weightsOnEdge_lsqr, & !$acc kiteAreasOnVertex, & !$acc edgeTangentVectors, & !$acc edgeNormalVectors, & diff --git a/components/mpas-ocean/src/shared/mpas_ocn_vel_hmix.F b/components/mpas-ocean/src/shared/mpas_ocn_vel_hmix.F index 0333fe66809a..86abc6a6554c 100644 --- a/components/mpas-ocean/src/shared/mpas_ocn_vel_hmix.F +++ b/components/mpas-ocean/src/shared/mpas_ocn_vel_hmix.F @@ -125,13 +125,16 @@ subroutine ocn_vel_hmix_tend(div, relVort, tend, err)!{{{ ! note that the user can choose multiple options and the ! tendencies will be added together + call ocn_vel_hmix_aust_del2_tend(div, tend, err1) call ocn_vel_hmix_del2_tend(div, relVort, tend, err1) err = ior(err1, err) call ocn_vel_hmix_leith_tend(div, relVort, tend, err1) err = ior(err1, err) + call ocn_vel_hmix_aust_del4_tend(div, tend, err1) call ocn_vel_hmix_del4_tend(div, relVort, tend, err1) + !call ocn_vel_hmix_div_del4_tend(div, tend, err1) err = ior(err1, err) call mpas_timer_stop("vel hmix") diff --git a/components/mpas-ocean/src/shared/mpas_ocn_vel_hmix_del2.F b/components/mpas-ocean/src/shared/mpas_ocn_vel_hmix_del2.F index 42fad6a10713..1dfa6888c4b7 100644 --- a/components/mpas-ocean/src/shared/mpas_ocn_vel_hmix_del2.F +++ b/components/mpas-ocean/src/shared/mpas_ocn_vel_hmix_del2.F @@ -42,6 +42,7 @@ module ocn_vel_hmix_del2 !-------------------------------------------------------------------- public :: ocn_vel_hmix_del2_tend, & + ocn_vel_hmix_aust_del2_tend, & ocn_vel_hmix_del2_init !------------------------------------------------------------------- @@ -51,10 +52,12 @@ module ocn_vel_hmix_del2 !-------------------------------------------------------------------- logical :: & - hmixDel2Off !< on/off flag to determine whether del2 chosen + hmixDel2Off, & !< on/off flag to determine whether del2 chosen + hmixDel2AustOff !< on/off flag for AUST closure real (kind=RKIND) :: & - viscDel2 !< viscosity coefficient for del2 horz mixing + viscDel2, & !< viscosity coefficient for del2 horz mixing + viscDel2Aust !< viscosity coefficient for aust !*********************************************************************** @@ -179,6 +182,111 @@ subroutine ocn_vel_hmix_del2_tend(div, relVort, tend, err)!{{{ end subroutine ocn_vel_hmix_del2_tend!}}} +!*********************************************************************** +! +! routine ocn_vel_hmix_aust_del2_tend +! +!> \brief Computes tendency for Laplacian horizontal momentum mixing +!> \author Sara Calandrini +!> \date March 2022 +!> \details +!> This routine computes the horizontal mixing tendency for momentum +!> based on a Laplacian form for the mixing, \f$\nu_2 \nabla^2 u\f$ +!> This tendency takes the form +!> \f$\nu( \nabla divergence)\f$, where \f$\nu\f$ is a viscosity, so +!> the vorticity part has been eliminated compared to the above +!> subroutine. This form is strictly only valid for constant \f$\nu\f$ . +! +!----------------------------------------------------------------------- + + subroutine ocn_vel_hmix_aust_del2_tend(div, tend, err)!{{{ + + !----------------------------------------------------------------- + ! input variables + !----------------------------------------------------------------- + + real (kind=RKIND), dimension(:,:), intent(in) :: & + div !< [in] velocity divergence + + !----------------------------------------------------------------- + ! input /output variables + !----------------------------------------------------------------- + + real (kind=RKIND), dimension(:,:), intent(inout) :: & + tend !< [inout] accumulated velocity tendency + + !----------------------------------------------------------------- + ! output variables + !----------------------------------------------------------------- + + integer, intent(out) :: err !< [out] error flag + + !----------------------------------------------------------------- + ! local variables + !----------------------------------------------------------------- + + integer :: & + iEdge, k, &! loop indices for edge, vertical loops + cell1, cell2 ! neighbor cell addresses across edge + + real (kind=RKIND) :: & + uDiff, &! velocity diffusion + dcEdgeInv, &! 1/dcEdge + visc2 ! scaled viscosity coefficient + + !----------------------------------------------------------------- + ! + ! exit if this mixing is not selected + ! + !----------------------------------------------------------------- + + err = 0 + + if (hmixDel2AustOff) return + + call mpas_timer_start("vel div2") + +#ifdef MPAS_OPENACC + !$acc parallel loop & + !$acc present(cellsOnEdge, maxLevelEdgeTop, minLevelEdgeBot, & + !$acc edgeMask, dcEdge, div, & + !$acc meshScalingDel2, tend) & + !$acc private(k, cell1, cell2, uDiff, & + !$acc dcEdgeInv, visc2) +#else + !$omp parallel + !$omp do schedule(runtime) & + !$omp private(k, cell1, cell2, uDiff, & + !$omp dcEdgeInv, visc2) +#endif + do iEdge = 1, nEdgesOwned + cell1 = cellsOnEdge(1,iEdge) + cell2 = cellsOnEdge(2,iEdge) + + dcEdgeInv = 1.0_RKIND / dcEdge(iEdge) + + visc2 = viscDel2Aust*meshScalingDel2(iEdge) + + do k = minLevelEdgeBot(iEdge), maxLevelEdgeTop(iEdge) + + uDiff = (div(k,cell2) - div(k,cell1))*dcEdgeInv + + tend(k,iEdge) = tend(k,iEdge) + & + edgeMask(k,iEdge)*visc2*uDiff + + end do + end do +#ifndef MPAS_OPENACC + !$omp end do + !$omp end parallel +#endif + + call mpas_timer_stop("vel div2") + + !-------------------------------------------------------------------- + + end subroutine ocn_vel_hmix_aust_del2_tend!}}} + !*********************************************************************** ! ! routine ocn_vel_hmix_del2_init @@ -206,8 +314,10 @@ subroutine ocn_vel_hmix_del2_init(err)!{{{ !*** Initialize return code and default variables err = 0 - hmixDel2Off = .true. - viscDel2 = 0.0_RKIND + hmixDel2Off = .true. + hmixDel2AustOff = .true. + viscDel2 = 0.0_RKIND + viscDel2Aust = 0.0_RKIND !*** reset values based on input configuration @@ -221,8 +331,22 @@ subroutine ocn_vel_hmix_del2_init(err)!{{{ MPAS_LOG_ERR) err = -1 endif + endif + if (config_use_mom_div2) then + hmixDel2AustOff = .false. + viscDel2Aust = config_mom_div2 + + if (viscDel2Aust <= 0.0_RKIND) then + call mpas_log_write( & + 'vel_hmix_del2_init: viscosity must be > 0', & + MPAS_LOG_ERR) + err = -1 + endif + + end if + !----------------------------------------------------------------- end subroutine ocn_vel_hmix_del2_init!}}} diff --git a/components/mpas-ocean/src/shared/mpas_ocn_vel_hmix_del4.F b/components/mpas-ocean/src/shared/mpas_ocn_vel_hmix_del4.F index 11c922ce6078..c7c2b5267ead 100644 --- a/components/mpas-ocean/src/shared/mpas_ocn_vel_hmix_del4.F +++ b/components/mpas-ocean/src/shared/mpas_ocn_vel_hmix_del4.F @@ -43,6 +43,8 @@ module ocn_vel_hmix_del4 !-------------------------------------------------------------------- public :: ocn_vel_hmix_del4_tend, & + ocn_vel_hmix_div_del4_tend, & + ocn_vel_hmix_aust_del4_tend, & ocn_vel_hmix_del4_init !-------------------------------------------------------------------- @@ -52,10 +54,12 @@ module ocn_vel_hmix_del4 !-------------------------------------------------------------------- logical :: & - hmixDel4Off !< on/off flag to determine whether del4 chosen + hmixDel4Off, & !< on/off flag to determine whether del4 chosen + hmixDel4AustOff !< on/off flag for AUST closure real (kind=RKIND) :: & viscDel4, &!< biharmonic viscosity coefficient + viscDel4Aust, &!< AUST viscosity coefficient divFactor !< factor for divergence term !*********************************************************************** @@ -307,6 +311,394 @@ subroutine ocn_vel_hmix_del4_tend(div, relVort, tend, err)!{{{ end subroutine ocn_vel_hmix_del4_tend!}}} +!*********************************************************************** +! +! routine ocn_vel_hmix_div_del4_tend +! +!> \brief Computes tendency for biharmonic horizontal momentum mixing +!> \author Sara Calandrini +!> \date March 2022 +!> \details +!> This routine computes the horizontal mixing tendency for momentum +!> based on a biharmonic form for the mixing. The del4 operator is +!> applied on the divergence part only, so the difference with the +!> previous subroutine is that any contribution from the vorticity +!> is eliminated. +! +!----------------------------------------------------------------------- + + subroutine ocn_vel_hmix_div_del4_tend(div, tend, err)!{{{ + + !----------------------------------------------------------------- + ! input variables + !----------------------------------------------------------------- + + real (kind=RKIND), dimension(:,:), intent(in) :: & + div !< [in] velocity divergence + + !----------------------------------------------------------------- + ! input/output variables + !----------------------------------------------------------------- + + real (kind=RKIND), dimension(:,:), intent(inout) :: & + tend !< [inout] accumulated velocity tendency + + !----------------------------------------------------------------- + ! output variables + !----------------------------------------------------------------- + + integer, intent(out) :: err !< [out] error flag + + !----------------------------------------------------------------- + ! local variables + !----------------------------------------------------------------- + + integer :: & + iEdge, iCell, iVertex, k, i, &! loop counters + cell1, cell2, &! neighbor cell addresses across edge + nEdges, nCells, nVertices ! num edges,cells,vertices incl halos + + real (kind=RKIND) :: & + uDiff, &! diffusion operator temporary + areaCellInv, &! 1/area of cell + areaTriInv, &! 1/area of triangle + dcEdgeInv, &! 1/dcEdge + visc4 ! scaled biharmonic viscosity coeff + + ! Scratch Arrays + real (kind=RKIND), dimension(:,:), allocatable :: & + del2div, &! + del2u + + integer :: kmin, kmax + + ! End preamble + !----------------------------------------------------------------- + ! Begin code + + !*** initialize error code and return if not selected + !*** otherwise start timer + + err = 0 + if (hmixDel4Off) return + call mpas_timer_start("vel del4") + + !*** allocate temporaries + + allocate(del2u(nVertLevels, nEdgesAll + 1), & + del2div(nVertLevels, nCellsAll + 1)) + !$acc enter data create(del2u, del2div) + + nEdges = nEdgesHalo(2) + + !Compute Laplacian of velocity +#ifdef MPAS_OPENACC + !$acc parallel loop & + !$acc present(cellsOnEdge, verticesOnEdge, minLevelEdgeBot, & + !$acc maxLevelEdgeTop, dcEdge, del2u, div) & + !$acc private(k, cell1, cell2, & + !$acc dcEdgeInv, kmin, kmax) +#else + !$omp parallel + !$omp do schedule(runtime) & + !$omp private(k, cell1, cell2, & + !$omp dcEdgeInv, kmin, kmax) +#endif + do iEdge = 1, nEdges + del2u(:, iEdge) = 0.0_RKIND + cell1 = cellsOnEdge(1,iEdge) + cell2 = cellsOnEdge(2,iEdge) + + kmin = minLevelEdgeBot(iEdge) + kmax = maxLevelEdgeTop(iEdge) + + dcEdgeInv = 1.0_RKIND / dcEdge(iEdge) + + do k=kmin, kmax + ! Compute \nabla^2 u = \nabla divergence + del2u(k,iEdge) = (div(k,cell2) - div(k,cell1))*dcEdgeInv + end do + end do +#ifndef MPAS_OPENACC + !$omp end do + !$omp end parallel +#endif + + nCells = nCellsHalo(1) + + ! Compute del2 of divergence +#ifdef MPAS_OPENACC + !$acc parallel loop & + !$acc present(nEdgesOnCell, edgesOnCell, minLevelCell, maxLevelCell, & + !$acc dvEdge, edgeSignOnCell, areaCell, del2u, del2div) & + !$acc private(i, k, iEdge, areaCellInv, kmin, kmax) +#else + !$omp parallel + !$omp do schedule(runtime) private(i, k, iEdge, areaCellInv, kmin, kmax) +#endif + do iCell = 1, nCells + kmin = minLevelCell(iCell) + kmax = maxLevelCell(iCell) + del2div(:, iCell) = 0.0_RKIND + areaCellInv = 1.0_RKIND / areaCell(iCell) + do i = 1, nEdgesOnCell(iCell) + iEdge = edgesOnCell(i, iCell) + do k = kmin, kmax + del2div(k,iCell) = del2div(k,iCell) - & + edgeSignOnCell(i,iCell)*dvEdge(iEdge) & + *del2u(k,iEdge)*areaCellInv + end do + end do + end do +#ifndef MPAS_OPENACC + !$omp end do + !$omp end parallel +#endif + + ! Compute final tendency \kappa \nabla^4 u + ! as \nabla div(\nabla^2 u) + +#ifdef MPAS_OPENACC + !$acc parallel loop & + !$acc present(cellsOnEdge, minLevelEdgeBot, maxLevelEdgeTop, & + !$acc dcEdge, meshScalingDel4, edgeMask, & + !$acc del2div, tend) & + !$acc private(k, cell1, cell2, dcEdgeInv, & + !$acc visc4, uDiff, kmin, kmax) +#else + !$omp parallel + !$omp do schedule(runtime) & + !$omp private(k, cell1, cell2, dcEdgeInv, & + !$omp visc4, uDiff, kmin, kmax) +#endif + do iEdge = 1, nEdgesOwned + kmin = minLevelEdgeBot(iEdge) + kmax = maxLevelEdgeTop(iEdge) + + cell1 = cellsOnEdge(1,iEdge) + cell2 = cellsOnEdge(2,iEdge) + + dcEdgeInv = 1.0_RKIND / dcEdge(iEdge) + visc4 = viscDel4*meshScalingDel4(iEdge) + + do k = kmin, kmax + uDiff = (del2div(k,cell2) - del2div(k,cell1))* & + dcEdgeInv + + tend(k,iEdge) = tend(k,iEdge) - & + edgeMask(k,iEdge)*uDiff*visc4 + + end do + end do +#ifndef MPAS_OPENACC + !$omp end do + !$omp end parallel +#endif + + !$acc exit data delete(del2u, del2div) + deallocate(del2u, & + del2div) + + call mpas_timer_stop("vel del4") + + !-------------------------------------------------------------------- + + end subroutine ocn_vel_hmix_div_del4_tend!}}} + +!*********************************************************************** +! +! routine ocn_vel_hmix_aust_del4_tend +! +!> \brief Computes tendency for biharmonic horizontal momentum mixing +!> \author Sara Calandrini +!> \date March 2022 +!> \details +!> This routine computes the horizontal mixing tendency for momentum +!> based on a biharmonic form for the mixing. The del4 operator is +!> applied on the divergence part only, and the visc4 is inserted +!> within the gradient operator. +! +!----------------------------------------------------------------------- + + subroutine ocn_vel_hmix_aust_del4_tend(div, tend, err)!{{{ + + !----------------------------------------------------------------- + ! input variables + !----------------------------------------------------------------- + + real (kind=RKIND), dimension(:,:), intent(in) :: & + div !< [in] velocity divergence + + !----------------------------------------------------------------- + ! input/output variables + !----------------------------------------------------------------- + + real (kind=RKIND), dimension(:,:), intent(inout) :: & + tend !< [inout] accumulated velocity tendency + + !----------------------------------------------------------------- + ! output variables + !----------------------------------------------------------------- + + integer, intent(out) :: err !< [out] error flag + + !----------------------------------------------------------------- + ! local variables + !----------------------------------------------------------------- + + integer :: & + iEdge, iCell, iVertex, k, i, &! loop counters + cell1, cell2, &! neighbor cell addresses across edge + nEdges, nCells, nVertices ! num edges,cells,vertices incl halos + + real (kind=RKIND) :: & + uDiff, &! diffusion operator temporary + areaCellInv, &! 1/area of cell + areaTriInv, &! 1/area of triangle + dcEdgeInv, &! 1/dcEdge + visc4 ! scaled biharmonic viscosity coeff + + ! Scratch Arrays + real (kind=RKIND), dimension(:,:), allocatable :: & + del2div, &! + del2u + + integer :: kmin, kmax + + ! End preamble + !----------------------------------------------------------------- + ! Begin code + + !*** initialize error code and return if not selected + !*** otherwise start timer + + err = 0 + if (hmixDel4AustOff) return + call mpas_timer_start("vel div4") + + !*** allocate temporaries + + allocate(del2u(nVertLevels, nEdgesAll + 1), & + del2div(nVertLevels, nCellsAll + 1)) + !$acc enter data create(del2u, del2div) + + nEdges = nEdgesHalo(2) + + !Compute Laplacian of velocity +#ifdef MPAS_OPENACC + !$acc parallel loop & + !$acc present(cellsOnEdge, verticesOnEdge, minLevelEdgeBot, & + !$acc maxLevelEdgeTop, dcEdge, del2u, div) & + !$acc private(k, cell1, cell2, visc4, & + !$acc dcEdgeInv, kmin, kmax) +#else + !$omp parallel + !$omp do schedule(runtime) & + !$omp private(k, cell1, cell2, visc4, & + !$omp dcEdgeInv, kmin, kmax) +#endif + do iEdge = 1, nEdges + del2u(:, iEdge) = 0.0_RKIND + cell1 = cellsOnEdge(1,iEdge) + cell2 = cellsOnEdge(2,iEdge) + + kmin = minLevelEdgeBot(iEdge) + kmax = maxLevelEdgeTop(iEdge) + + dcEdgeInv = 1.0_RKIND / dcEdge(iEdge) + visc4 = viscDel4Aust*meshScalingDel4(iEdge) + + do k=kmin, kmax + ! Compute \nabla^2 u = \nabla sqrt(visc4)*divergence + del2u(k,iEdge) = sqrt(visc4)*(div(k,cell2) - div(k,cell1))*dcEdgeInv + end do + end do +#ifndef MPAS_OPENACC + !$omp end do + !$omp end parallel +#endif + + nCells = nCellsHalo(1) + + ! Compute del2 of divergence +#ifdef MPAS_OPENACC + !$acc parallel loop & + !$acc present(nEdgesOnCell, edgesOnCell, minLevelCell, maxLevelCell, & + !$acc dvEdge, edgeSignOnCell, areaCell, del2u, del2div) & + !$acc private(i, k, iEdge, areaCellInv, kmin, kmax) +#else + !$omp parallel + !$omp do schedule(runtime) private(i, k, iEdge, areaCellInv, kmin, kmax) +#endif + do iCell = 1, nCells + kmin = minLevelCell(iCell) + kmax = maxLevelCell(iCell) + del2div(:, iCell) = 0.0_RKIND + areaCellInv = 1.0_RKIND / areaCell(iCell) + do i = 1, nEdgesOnCell(iCell) + iEdge = edgesOnCell(i, iCell) + do k = kmin, kmax + del2div(k,iCell) = del2div(k,iCell) - & + edgeSignOnCell(i,iCell)*dvEdge(iEdge) & + *del2u(k,iEdge)*areaCellInv + end do + end do + end do +#ifndef MPAS_OPENACC + !$omp end do + !$omp end parallel +#endif + + ! Compute final tendency \kappa \nabla^4 u + ! as \nabla div(\nabla^2 u) + +#ifdef MPAS_OPENACC + !$acc parallel loop & + !$acc present(cellsOnEdge, minLevelEdgeBot, maxLevelEdgeTop, & + !$acc dcEdge, meshScalingDel4, edgeMask, & + !$acc del2div, tend) & + !$acc private(k, cell1, cell2, dcEdgeInv, & + !$acc visc4, uDiff, kmin, kmax) +#else + !$omp parallel + !$omp do schedule(runtime) & + !$omp private(k, cell1, cell2, dcEdgeInv, & + !$omp visc4, uDiff, kmin, kmax) +#endif + do iEdge = 1, nEdgesOwned + kmin = minLevelEdgeBot(iEdge) + kmax = maxLevelEdgeTop(iEdge) + + cell1 = cellsOnEdge(1,iEdge) + cell2 = cellsOnEdge(2,iEdge) + + dcEdgeInv = 1.0_RKIND / dcEdge(iEdge) + visc4 = viscDel4Aust*meshScalingDel4(iEdge) + + do k = kmin, kmax + uDiff = sqrt(visc4)*(del2div(k,cell2) - del2div(k,cell1))* & + dcEdgeInv + + tend(k,iEdge) = tend(k,iEdge) - & + edgeMask(k,iEdge)*uDiff + + end do + end do +#ifndef MPAS_OPENACC + !$omp end do + !$omp end parallel +#endif + + !$acc exit data delete(del2u, del2div) + deallocate(del2u, & + del2div) + + call mpas_timer_stop("vel div4") + + !-------------------------------------------------------------------- + + end subroutine ocn_vel_hmix_aust_del4_tend!}}} + !*********************************************************************** ! ! routine ocn_vel_hmix_del4_init @@ -335,16 +727,18 @@ subroutine ocn_vel_hmix_del4_init(err)!{{{ !*** initialize return error code and set module variable defaults err = 0 - hmixDel4Off = .true. - viscDel4 = 0.0_RKIND - divFactor = 1.0_RKIND + hmixDel4Off = .true. + hmixDel4AustOff = .true. + viscDel4 = 0.0_RKIND + viscDel4Aust = 0.0_RKIND + divFactor = 1.0_RKIND !*** reset module variables based on input configuration if (config_use_mom_del4) then - hmixDel4Off = .false. - viscDel4 = config_mom_del4 - divFactor = config_mom_del4_div_factor + hmixDel4Off = .false. + viscDel4 = config_mom_del4 + divFactor = config_mom_del4_div_factor if ( config_mom_del4 <= 0.0_RKIND ) then call mpas_log_write( & 'vel_hmix_del4_init: config_mom_del4 must be > 0 ', & @@ -353,6 +747,17 @@ subroutine ocn_vel_hmix_del4_init(err)!{{{ endif endif + if (config_use_mom_div4) then + hmixDel4AustOff = .false. + viscDel4Aust = config_mom_div4 + if ( config_mom_div4 <= 0.0_RKIND ) then + call mpas_log_write( & + 'vel_hmix_del4_init: config_mom_div4 must be > 0 ', & + MPAS_LOG_ERR) + err = -1 + endif + endif + !-------------------------------------------------------------------- end subroutine ocn_vel_hmix_del4_init!}}}