<p><b>dwj07@fsu.edu</b> 2011-09-15 15:59:53 -0600 (Thu, 15 Sep 2011)</p><p><br>
-- Branch commit --<br>
<br>
Adding modules for Coriolis force, Pressure Gradient, Vertical advection, and Horizontal mixing on the momentum equation.<br>
</p><hr noshade><pre><font color="gray">Modified: branches/ocean_projects/performance/src/core_ocean/Makefile
===================================================================
--- branches/ocean_projects/performance/src/core_ocean/Makefile 2011-09-15 17:43:05 UTC (rev 1001)
+++ branches/ocean_projects/performance/src/core_ocean/Makefile 2011-09-15 21:59:53 UTC (rev 1002)
@@ -3,6 +3,12 @@
OBJS = module_mpas_core.o \
module_test_cases.o \
module_advection.o \
+ module_OcnVelCoriolis.o \
+ module_OcnVelVadv.o \
+ module_OcnVelHmix.o \
+ module_OcnVelHmixDel2.o \
+ module_OcnVelHmixDel4.o \
+ module_OcnVelPressureGrad.o \
module_time_integration.o \
module_global_diagnostics.o
@@ -19,8 +25,29 @@
module_global_diagnostics.o:
-module_mpas_core.o: module_advection.o module_global_diagnostics.o module_test_cases.o module_time_integration.o
+module_OcnVelPressureGrad.o:
+module_OcnVelVadv.o:
+
+module_OcnVelHmix.o: module_OcnVelHmixDel2.o module_OcnVelHmixDel4.o
+
+module_OcnVelHmixDel2.o:
+
+module_OcnVelHmixDel4.o:
+
+module_OcnVelCoriolis.o:
+
+module_mpas_core.o: module_advection.o \
+ module_global_diagnostics.o \
+ module_test_cases.o \
+ module_OcnVelCoriolis.o \
+ module_OcnVelVadv.o \
+ module_OcnVelHmix.o \
+ module_OcnVelHmixDel2.o \
+ module_OcnVelHmixDel4.o \
+ module_OcnVelPressureGrad.o \
+ module_time_integration.o
+
clean:
$(RM) *.o *.mod *.f90 libdycore.a
Added: branches/ocean_projects/performance/src/core_ocean/module_OcnVelCoriolis.F
===================================================================
--- branches/ocean_projects/performance/src/core_ocean/module_OcnVelCoriolis.F (rev 0)
+++ branches/ocean_projects/performance/src/core_ocean/module_OcnVelCoriolis.F 2011-09-15 21:59:53 UTC (rev 1002)
@@ -0,0 +1,186 @@
+!|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
+!
+! OcnVelCoriolis
+!
+!> \brief MPAS ocean horizontal momentum mixing driver
+!> \author Doug Jacobsen
+!> \date 15 September 2011
+!> \version SVN:$Id:$
+!> \details
+!> This module contains the routine for computing
+!> tendencies from the coriolis force.
+!>
+!
+!-----------------------------------------------------------------------
+
+module OcnVelCoriolis
+
+ use grid_types
+ use configure
+
+ implicit none
+ private
+ save
+
+ !--------------------------------------------------------------------
+ !
+ ! Public parameters
+ !
+ !--------------------------------------------------------------------
+
+ !--------------------------------------------------------------------
+ !
+ ! Public member functions
+ !
+ !--------------------------------------------------------------------
+
+ public :: OcnVelCoriolisTend, &
+ OcnVelCoriolisInit
+
+ !--------------------------------------------------------------------
+ !
+ ! Private module variables
+ !
+ !--------------------------------------------------------------------
+
+
+!***********************************************************************
+
+contains
+
+!***********************************************************************
+!
+! routine OcnVelCoriolisTend
+!
+!> \brief Computes tendency term for coriolis force
+!> \author Doug Jacobsen
+!> \date 15 September 2011
+!> \version SVN:$Id$
+!> \details
+!> This routine computes the coriolis tendency for momentum
+!> based on current state.
+!
+!-----------------------------------------------------------------------
+
+ subroutine OcnVelCoriolisTend(grid, pv_edge, h_edge, u, ke, tend, err)
+
+ !-----------------------------------------------------------------
+ !
+ ! input variables
+ !
+ !-----------------------------------------------------------------
+
+ real (kind=RKIND), dimension(:,:), intent(in) :: &
+ pv_edge, h_edge, u, ke
+
+ type (mesh_type), intent(in) :: &
+ grid !< Input: grid information
+
+ !-----------------------------------------------------------------
+ !
+ ! input/output variables
+ !
+ !-----------------------------------------------------------------
+
+ real (kind=RKIND), dimension(:,:), intent(inout) :: &
+ tend !< Input/Output: velocity tendency
+
+ !-----------------------------------------------------------------
+ !
+ ! output variables
+ !
+ !-----------------------------------------------------------------
+
+ integer, intent(out) :: err !< Output: error flag
+
+ !-----------------------------------------------------------------
+ !
+ ! local variables
+ !
+ !-----------------------------------------------------------------
+
+ integer, dimension(:), pointer :: maxLevelEdgeTop, nEdgesOnEdge
+ integer, dimension(:,:), pointer :: cellsOnEdge, edgesOnEdge
+ real (kind=RKIND), dimension(:,:), pointer :: weightsOnEdge
+ real (kind=RKIND), dimension(:), pointer :: dcEdge
+
+ integer :: j, k
+ integer :: cell1, cell2, nEdgesSolve, iEdge, eoe
+ real (kind=RKIND) :: workpv, q
+
+ err = 0
+
+ maxLevelEdgeTop => grid % maxLevelEdgeTop % array
+ nEdgesOnEdge => grid % nEdgesOnEdge % array
+ cellsOnEdge => grid % cellsOnEdge % array
+ edgesOnEdge => grid % edgesOnEdge % array
+ weightsOnEdge => grid % weightsOnEdge % array
+ dcEdge => grid % dcEdge % array
+
+ nEdgesSolve = grid % nEdgesSolve
+
+ do iEdge=1,nEdgesSolve
+ cell1 = cellsOnEdge(1,iEdge)
+ cell2 = cellsOnEdge(2,iEdge)
+
+ do k=1,maxLevelEdgeTop(iEdge)
+
+ q = 0.0
+ do j = 1,nEdgesOnEdge(iEdge)
+ eoe = edgesOnEdge(j,iEdge)
+ workpv = 0.5 * (pv_edge(k,iEdge) + pv_edge(k,eoe))
+ q = q + weightsOnEdge(j,iEdge) * u(k,eoe) * workpv * h_edge(k,eoe)
+ end do
+
+ tend(k,iEdge) = tend(k,iEdge) &
+ + q &
+ - ( ke(k,cell2) - ke(k,cell1) ) / dcEdge(iEdge)
+
+ end do
+ end do
+
+
+ !--------------------------------------------------------------------
+
+ end subroutine OcnVelCoriolisTend
+
+!***********************************************************************
+!
+! routine OcnVelCoriolisInit
+!
+!> \brief Initializes ocean momentum horizontal mixing quantities
+!> \author Doug Jacobsen
+!> \date 15 September 2011
+!> \version SVN:$Id$
+!> \details
+!> This routine initializes a variety of quantities related to
+!> horizontal velocity mixing in the ocean. Since a variety of
+!> parameterizations are available, this routine primarily calls the
+!> individual init routines for each parameterization.
+!
+!-----------------------------------------------------------------------
+
+
+ subroutine OcnVelCoriolisInit(err)
+
+ !--------------------------------------------------------------------
+
+ !-----------------------------------------------------------------
+ !
+ ! Output variables
+ !
+ !-----------------------------------------------------------------
+
+ integer, intent(out) :: err
+
+ err = 0
+
+ !--------------------------------------------------------------------
+
+ end subroutine OcnVelCoriolisInit
+
+!***********************************************************************
+
+end module OcnVelCoriolis
+
+!|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
Added: branches/ocean_projects/performance/src/core_ocean/module_OcnVelHmix.F
===================================================================
--- branches/ocean_projects/performance/src/core_ocean/module_OcnVelHmix.F (rev 0)
+++ branches/ocean_projects/performance/src/core_ocean/module_OcnVelHmix.F 2011-09-15 21:59:53 UTC (rev 1002)
@@ -0,0 +1,175 @@
+!|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
+!
+! OcnVelHmix
+!
+!> \brief MPAS ocean horizontal momentum mixing driver
+!> \author Phil Jones, Doug Jacobsen
+!> \date 15 September 2011
+!> \version SVN:$Id:$
+!> \details
+!> This module contains the main driver routine for computing
+!> horizontal mixing tendencies.
+!>
+!> It provides an init and a tend function. Each are described below.
+!
+!-----------------------------------------------------------------------
+
+module OcnVelHmix
+
+ use grid_types
+ use configure
+ use OcnVelHmixDel2
+ use OcnVelHmixDel4
+
+ implicit none
+ private
+ save
+
+ !--------------------------------------------------------------------
+ !
+ ! Public parameters
+ !
+ !--------------------------------------------------------------------
+
+ !--------------------------------------------------------------------
+ !
+ ! Public member functions
+ !
+ !--------------------------------------------------------------------
+
+ public :: OcnVelHmixTend, &
+ OcnVelHmixInit
+
+ !--------------------------------------------------------------------
+ !
+ ! Private module variables
+ !
+ !--------------------------------------------------------------------
+
+
+!***********************************************************************
+
+contains
+
+!***********************************************************************
+!
+! routine OcnVelHmixTend
+!
+!> \brief Computes tendency term for horizontal momentum mixing
+!> \author Phil Jones, Doug Jacobsen
+!> \date 15 September 2011
+!> \version SVN:$Id$
+!> \details
+!> This routine computes the horizontal mixing tendency for momentum
+!> based on current state and user choices of mixing parameterization.
+!> Multiple parameterizations may be chosen and added together. These
+!> tendencies are generally computed by calling the specific routine
+!> for the chosen parameterization, so this routine is primarily a
+!> driver for managing these choices.
+!
+!-----------------------------------------------------------------------
+
+ subroutine OcnVelHmixTend(grid, divergence, vorticity, tend, err)
+
+ !-----------------------------------------------------------------
+ !
+ ! input variables
+ !
+ !-----------------------------------------------------------------
+
+ real (kind=RKIND), dimension(:,:), intent(in) :: &
+ divergence !< Input: velocity divergence
+
+ real (kind=RKIND), dimension(:,:), intent(in) :: &
+ vorticity !< Input: vorticity
+
+ type (mesh_type), intent(in) :: &
+ grid !< Input: grid information
+
+ !-----------------------------------------------------------------
+ !
+ ! input/output variables
+ !
+ !-----------------------------------------------------------------
+
+ real (kind=RKIND), dimension(:,:), intent(inout) :: &
+ tend !< Input/Output: velocity tendency
+
+ !-----------------------------------------------------------------
+ !
+ ! output variables
+ !
+ !-----------------------------------------------------------------
+
+ integer, intent(out) :: err
+
+ !-----------------------------------------------------------------
+ !
+ ! local variables
+ !
+ !-----------------------------------------------------------------
+
+ integer :: err1, err2
+
+ !-----------------------------------------------------------------
+ !
+ ! call relevant routines for computing tendencies
+ ! note that the user can choose multiple options and the
+ ! tendencies will be added together
+ !
+ !-----------------------------------------------------------------
+
+ call OcnVelHmixDel2Tend(grid, divergence, vorticity, tend, err1)
+ call OcnVelHmixDel4Tend(grid, divergence, vorticity, tend, err2)
+
+ err = err1 .or. err2
+
+ !--------------------------------------------------------------------
+
+ end subroutine OcnVelHmixTend
+
+!***********************************************************************
+!
+! routine OcnVelHmixInit
+!
+!> \brief Initializes ocean momentum horizontal mixing quantities
+!> \author Phil Jones, Doug Jacobsen
+!> \date 15 September 2011
+!> \version SVN:$Id$
+!> \details
+!> This routine initializes a variety of quantities related to
+!> horizontal velocity mixing in the ocean. Since a variety of
+!> parameterizations are available, this routine primarily calls the
+!> individual init routines for each parameterization.
+!
+!-----------------------------------------------------------------------
+
+
+ subroutine OcnVelHmixInit(err)
+
+ !--------------------------------------------------------------------
+
+ !-----------------------------------------------------------------
+ !
+ ! call individual init routines for each parameterization
+ !
+ !-----------------------------------------------------------------
+
+ integer, intent(out) :: err
+
+ integer :: err1, err2
+
+ call OcnVelHmixDel2Init(err1)
+ call OcnVelHmixDel4Init(err2)
+
+ err = err1 .or. err2
+
+ !--------------------------------------------------------------------
+
+ end subroutine OcnVelHmixInit
+
+!***********************************************************************
+
+end module OcnVelHmix
+
+!|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
Added: branches/ocean_projects/performance/src/core_ocean/module_OcnVelHmixDel2.F
===================================================================
--- branches/ocean_projects/performance/src/core_ocean/module_OcnVelHmixDel2.F (rev 0)
+++ branches/ocean_projects/performance/src/core_ocean/module_OcnVelHmixDel2.F 2011-09-15 21:59:53 UTC (rev 1002)
@@ -0,0 +1,224 @@
+!|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
+!
+! OcnVelHmixDel2
+!
+!> \brief Ocean horizontal mixing - Laplacian parameterization
+!> \author Phil Jones, Doug Jacobsen
+!> \date 15 September 2011
+!> \version SVN:$Id:$
+!> \details
+!> This module contains routines for computing horizontal mixing
+!> tendencies using a Laplacian formulation.
+!
+!-----------------------------------------------------------------------
+
+module OcnVelHmixDel2
+
+ use grid_types
+ use configure
+ use timer
+
+ implicit none
+ private
+ save
+
+ !--------------------------------------------------------------------
+ !
+ ! Public parameters
+ !
+ !--------------------------------------------------------------------
+
+ !--------------------------------------------------------------------
+ !
+ ! Public member functions
+ !
+ !--------------------------------------------------------------------
+
+ public :: OcnVelHmixDel2Tend, &
+ OcnVelHmixDel2Init
+
+ !-------------------------------------------------------------------
+ !
+ ! Private module variables
+ !
+ !--------------------------------------------------------------------
+
+ logical :: &
+ hmixDel2On !< local flag to determine whether del2 chosen
+
+ real (kind=RKIND) :: &
+ eddyVisc2, &!< base eddy diffusivity for Laplacian
+ viscVortCoef
+
+!***********************************************************************
+
+contains
+
+!***********************************************************************
+!
+! routine OcnVelHmixDel2Tend
+!
+!> \brief Computes tendency term for Laplacian horizontal momentum mixing
+!> \author Phil Jones, Doug Jacobsen
+!> \date 22 August 2011
+!> \version SVN:$Id$
+!> \details
+!> This routine computes the horizontal mixing tendency for momentum
+!> based on a Laplacian form for the mixing, </font>
<font color="black">u_2 </font>
<font color="blue">abla^2 u
+!> This tendency takes the
+!> form </font>
<font color="black">u( </font>
<font color="black">abla divergence + k \times </font>
<font color="blue">abla vorticity ),
+!> where </font>
<font color="blue">u is a viscosity and k is the vertical unit vector.
+!> This form is strictly only valid for constant </font>
<font color="blue">u .
+!
+!-----------------------------------------------------------------------
+
+ subroutine OcnVelHmixDel2Tend(grid, divergence, vorticity, tend, err)
+
+ !-----------------------------------------------------------------
+ !
+ ! input variables
+ !
+ !-----------------------------------------------------------------
+
+ real (kind=RKIND), dimension(:,:), intent(in) :: &
+ divergence !< Input: velocity divergence
+
+ real (kind=RKIND), dimension(:,:), intent(in) :: &
+ vorticity !< Input: vorticity
+
+ type (mesh_type), intent(in) :: &
+ grid !< Input: grid information
+
+ !-----------------------------------------------------------------
+ !
+ ! input/output variables
+ !
+ !-----------------------------------------------------------------
+
+ real (kind=RKIND), dimension(:,:), intent(inout) :: &
+ tend !< Input/Output: velocity tendency
+
+
+ !-----------------------------------------------------------------
+ !
+ ! output variables
+ !
+ !-----------------------------------------------------------------
+
+ integer, intent(out) :: err
+
+ !-----------------------------------------------------------------
+ !
+ ! local variables
+ !
+ !-----------------------------------------------------------------
+
+ integer :: iEdge, nEdgesSolve, cell1, cell2, vertex1, vertex2
+ integer :: k
+ integer, dimension(:), pointer :: maxLevelEdgeTop
+ integer, dimension(:,:), pointer :: cellsOnEdge, verticesOnEdge
+
+ real (kind=RKIND) :: u_diffusion
+ real (kind=RKIND), dimension(:), pointer :: meshScalingDel2, &
+ dcEdge, dvEdge
+
+ !-----------------------------------------------------------------
+ !
+ ! exit if this mixing is not selected
+ !
+ !-----------------------------------------------------------------
+
+ err = 0
+
+ if(.not.hmixDel2On) return
+
+ call timer_start("compute_tend_u-horiz mix-del2")
+
+ nEdgesSolve = grid % nEdgesSolve
+ maxLevelEdgeTop => grid % maxLevelEdgeTop % array
+ cellsOnEdge => grid % cellsOnEdge % array
+ verticesOnEdge => grid % verticesOnEdge % array
+ meshScalingDel2 => grid % meshScalingDel2 % array
+ dcEdge => grid % dcEdge % array
+ dvEdge => grid % dvEdge % array
+
+ do iEdge=1,nEdgesSolve
+ cell1 = cellsOnEdge(1,iEdge)
+ cell2 = cellsOnEdge(2,iEdge)
+ vertex1 = verticesOnEdge(1,iEdge)
+ vertex2 = verticesOnEdge(2,iEdge)
+
+ do k=1,maxLevelEdgeTop(iEdge)
+
+ ! Here -( vorticity(k,vertex2) - vorticity(k,vertex1) ) / dvEdge(iEdge)
+ ! is - </font>
<font color="blue">abla vorticity pointing from vertex 2 to vertex 1, or equivalently
+ ! + k \times </font>
<font color="gray">abla vorticity pointing from cell1 to cell2.
+
+ u_diffusion = ( divergence(k,cell2) - divergence(k,cell1) ) / dcEdge(iEdge) &
+ -viscVortCoef &
+ *( vorticity(k,vertex2) - vorticity(k,vertex1) ) / dvEdge(iEdge)
+
+ u_diffusion = meshScalingDel2(iEdge) * eddyVisc2 * u_diffusion
+
+ tend(k,iEdge) = tend(k,iEdge) + u_diffusion
+
+ end do
+ end do
+
+ call timer_stop("compute_tend_u-horiz mix-del2")
+
+ !--------------------------------------------------------------------
+
+ end subroutine OcnVelHmixDel2Tend
+
+!***********************************************************************
+!
+! routine OcnVelHmixDel2Init
+!
+!> \brief Initializes ocean momentum Laplacian horizontal mixing
+!> \author Phil Jones, Doug Jacobsen
+!> \date 15 September 2011
+!> \version SVN:$Id$
+!> \details
+!> This routine initializes a variety of quantities related to
+!> Laplacian horizontal momentum mixing in the ocean.
+!
+!-----------------------------------------------------------------------
+
+ subroutine OcnVelHmixDel2Init(err)
+
+
+ integer, intent(out) :: err
+
+ !--------------------------------------------------------------------
+ !
+ ! set some local module variables based on input config choices
+ !
+ !--------------------------------------------------------------------
+
+ err = 0
+
+ hmixDel2On = .false.
+
+ if ( config_h_mom_eddy_visc2 > 0.0 ) then
+ hmixDel2On = .true.
+ eddyVisc2 = config_h_mom_eddy_visc2
+
+
+ if (config_visc_vorticity_term) then
+ viscVortCoef = 1.0
+ else
+ viscVortCoef = 0.0
+ endif
+ endif
+
+
+ !--------------------------------------------------------------------
+
+ end subroutine OcnVelHmixDel2Init
+
+!***********************************************************************
+
+end module OcnVelHmixDel2
+
+!|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
Added: branches/ocean_projects/performance/src/core_ocean/module_OcnVelHmixDel4.F
===================================================================
--- branches/ocean_projects/performance/src/core_ocean/module_OcnVelHmixDel4.F (rev 0)
+++ branches/ocean_projects/performance/src/core_ocean/module_OcnVelHmixDel4.F 2011-09-15 21:59:53 UTC (rev 1002)
@@ -0,0 +1,300 @@
+!|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
+!
+! OcnVelHmixDel4
+!
+!> \brief Ocean horizontal mixing - biharmonic parameterization
+!> \author Phil Jones, Doug Jacobsen
+!> \date 15 September 2011
+!> \version SVN:$Id:$
+!> \details
+!> This module contains routines and variables for computing
+!> horizontal mixing tendencies using a biharmonic formulation.
+!
+!-----------------------------------------------------------------------
+
+module OcnVelHmixDel4
+
+ use grid_types
+ use configure
+ use timer
+
+ implicit none
+ private
+ save
+
+ !--------------------------------------------------------------------
+ !
+ ! Public parameters
+ !
+ !--------------------------------------------------------------------
+
+ !--------------------------------------------------------------------
+ !
+ ! Public member functions
+ !
+ !--------------------------------------------------------------------
+
+ public :: OcnVelHmixDel4Tend, &
+ OcnVelHmixDel4Init
+
+ !--------------------------------------------------------------------
+ !
+ ! Private module variables
+ !
+ !--------------------------------------------------------------------
+
+ logical :: &
+ hmixDel4On !< local flag to determine whether del4 chosen
+
+ real (kind=RKIND) :: &
+ eddyVisc4, &!< base eddy diffusivity for biharmonic
+ viscVortCoef
+
+
+!***********************************************************************
+
+contains
+
+!***********************************************************************
+!
+! routine OcnVelHmixDel4Tend
+!
+!> \brief Computes tendency term for biharmonic horizontal momentum mixing
+!> \author Phil Jones, Doug Jacobsen
+!> \date 15 September 2011
+!> \version SVN:$Id$
+!> \details
+!> This routine computes the horizontal mixing tendency for momentum
+!> based on a biharmonic form for the mixing. This mixing tendency
+!> takes the form -</font>
<font color="black">u_4 </font>
<font color="blue">abla^4 u
+!> but is computed as
+!> </font>
<font color="black">abla^2 u = </font>
<font color="black">abla divergence + k \times </font>
<font color="blue">abla vorticity
+!> applied recursively.
+!> This formulation is only valid for constant </font>
<font color="blue">u_4 .
+!
+!-----------------------------------------------------------------------
+
+ subroutine OcnVelHmixDel4Tend(grid, divergence, vorticity, tend, err)
+
+ !-----------------------------------------------------------------
+ !
+ ! input variables
+ !
+ !-----------------------------------------------------------------
+
+ real (kind=RKIND), dimension(:,:), intent(in) :: &
+ divergence !< Input: velocity divergence
+
+ real (kind=RKIND), dimension(:,:), intent(in) :: &
+ vorticity !< Input: vorticity
+
+ type (mesh_type), intent(in) :: &
+ grid !< Input: grid information
+
+ !-----------------------------------------------------------------
+ !
+ ! input/output variables
+ !
+ !-----------------------------------------------------------------
+
+ real (kind=RKIND), dimension(:,:), intent(inout) :: &
+ tend !< Input/Output: velocity tendency
+
+ !-----------------------------------------------------------------
+ !
+ ! output variables
+ !
+ !-----------------------------------------------------------------
+
+ integer, intent(out) :: err
+
+ !-----------------------------------------------------------------
+ !
+ ! local variables
+ !
+ !-----------------------------------------------------------------
+
+ integer :: iEdge, nEdges, cell1, cell2, vertex1, vertex2, k
+ integer :: iCell, iVertex
+ integer :: nVertices, nVertLevels, nCells
+
+ integer, dimension(:), pointer :: maxLevelEdgeTop, maxLevelVertexBot, &
+ maxLevelCell
+ integer, dimension(:,:), pointer :: cellsOnEdge, verticesOnEdge
+
+
+ real (kind=RKIND) :: u_diffusion, r
+ real (kind=RKIND), dimension(:), pointer :: dcEdge, dvEdge, areaTriangle, &
+ meshScalingDel4, areaCell
+
+ real (kind=RKIND), dimension(:,:), allocatable :: delsq_divergence, &
+ delsq_u, delsq_circulation, delsq_vorticity
+
+ err = 0
+
+ if(.not.hmixDel4On) return
+
+ call timer_start("compute_tend-horiz mix-del4")
+
+ nCells = grid % nCells
+ nEdges = grid % nEdges
+ nVertices = grid % nVertices
+ nVertLevels = grid % nVertLevels
+ maxLevelEdgeTop => grid % maxLevelEdgeTop % array
+ maxLevelVertexBot => grid % maxLevelVertexBot % array
+ maxLevelCell => grid % maxLevelCell % array
+ cellsOnEdge => grid % cellsOnEdge % array
+ verticesOnEdge => grid % verticesOnEdge % array
+ dcEdge => grid % dcEdge % array
+ dvEdge => grid % dvEdge % array
+ areaTriangle => grid % areaTriangle % array
+ areaCell => grid % areaCell % array
+ meshScalingDel4 => grid % meshScalingDel4 % array
+
+ allocate(delsq_divergence(nVertLevels, nCells+1))
+ allocate(delsq_u(nVertLevels, nEdges+1))
+ allocate(delsq_circulation(nVertLevels, nVertices+1))
+ allocate(delsq_vorticity(nVertLevels, nVertices+1))
+
+ delsq_u(:,:) = 0.0
+
+ ! Compute </font>
<font color="black">abla^2 u = </font>
<font color="black">abla divergence + k \times </font>
<font color="blue">abla vorticity
+ do iEdge=1,grid % nEdges
+ cell1 = cellsOnEdge(1,iEdge)
+ cell2 = cellsOnEdge(2,iEdge)
+ vertex1 = verticesOnEdge(1,iEdge)
+ vertex2 = verticesOnEdge(2,iEdge)
+
+ do k=1,maxLevelEdgeTop(iEdge)
+
+ delsq_u(k,iEdge) = &
+ ( divergence(k,cell2) - divergence(k,cell1) ) / dcEdge(iEdge) &
+ -viscVortCoef &
+ *( vorticity(k,vertex2) - vorticity(k,vertex1)) / dvEdge(iEdge)
+
+ end do
+ end do
+
+ ! vorticity using </font>
<font color="blue">abla^2 u
+ delsq_circulation(:,:) = 0.0
+ do iEdge=1,nEdges
+ vertex1 = verticesOnEdge(1,iEdge)
+ vertex2 = verticesOnEdge(2,iEdge)
+ do k=1,maxLevelEdgeTop(iEdge)
+ delsq_circulation(k,vertex1) = delsq_circulation(k,vertex1) &
+ - dcEdge(iEdge) * delsq_u(k,iEdge)
+ delsq_circulation(k,vertex2) = delsq_circulation(k,vertex2) &
+ + dcEdge(iEdge) * delsq_u(k,iEdge)
+ end do
+ end do
+ do iVertex=1,nVertices
+ r = 1.0 / areaTriangle(iVertex)
+ do k=1,maxLevelVertexBot(iVertex)
+ delsq_vorticity(k,iVertex) = delsq_circulation(k,iVertex) * r
+ end do
+ end do
+
+ ! Divergence using </font>
<font color="blue">abla^2 u
+ delsq_divergence(:,:) = 0.0
+ do iEdge=1,nEdges
+ cell1 = cellsOnEdge(1,iEdge)
+ cell2 = cellsOnEdge(2,iEdge)
+ do k=1,maxLevelEdgeTop(iEdge)
+ delsq_divergence(k,cell1) = delsq_divergence(k,cell1) &
+ + delsq_u(k,iEdge)*dvEdge(iEdge)
+ delsq_divergence(k,cell2) = delsq_divergence(k,cell2) &
+ - delsq_u(k,iEdge)*dvEdge(iEdge)
+ end do
+ end do
+ do iCell = 1,nCells
+ r = 1.0 / areaCell(iCell)
+ do k = 1,maxLevelCell(iCell)
+ delsq_divergence(k,iCell) = delsq_divergence(k,iCell) * r
+ end do
+ end do
+
+ ! Compute - \kappa </font>
<font color="blue">abla^4 u
+ ! as </font>
<font color="black">abla div(</font>
<font color="black">abla^2 u) + k \times </font>
<font color="black">abla ( k \cross curl(</font>
<font color="gray">abla^2 u) )
+ do iEdge=1,grid % nEdgesSolve
+ cell1 = cellsOnEdge(1,iEdge)
+ cell2 = cellsOnEdge(2,iEdge)
+ vertex1 = verticesOnEdge(1,iEdge)
+ vertex2 = verticesOnEdge(2,iEdge)
+
+ do k=1,maxLevelEdgeTop(iEdge)
+ delsq_u(k,iEdge) = &
+ ( divergence(k,cell2) - divergence(k,cell1) ) / dcEdge(iEdge) &
+ -( vorticity(k,vertex2) - vorticity(k,vertex1)) / dvEdge(iEdge)
+
+ u_diffusion = ( delsq_divergence(k,cell2) &
+ - delsq_divergence(k,cell1) ) / dcEdge(iEdge) &
+ -viscVortCoef &
+ *( delsq_vorticity(k,vertex2) &
+ - delsq_vorticity(k,vertex1) ) / dvEdge(iEdge)
+
+ u_diffusion = meshScalingDel4(iEdge) * eddyVisc4 * u_diffusion
+
+ tend(k,iEdge) = tend(k,iEdge) - u_diffusion
+ end do
+ end do
+
+ deallocate(delsq_divergence)
+ deallocate(delsq_u)
+ deallocate(delsq_circulation)
+ deallocate(delsq_vorticity)
+
+ call timer_stop("compute_tend-horiz mix-del4")
+
+ !--------------------------------------------------------------------
+
+ end subroutine OcnVelHmixDel4Tend
+
+!***********************************************************************
+!
+! routine OcnVelHmixDel4Init
+!
+!> \brief Initializes ocean momentum biharmonic horizontal mixing
+!> \author Phil Jones, Doug Jacobsen
+!> \date 15 September 2011
+!> \version SVN:$Id$
+!> \details
+!> This routine initializes a variety of quantities related to
+!> biharmonic horizontal tracer mixing in the ocean.
+!
+!-----------------------------------------------------------------------
+
+
+ subroutine OcnVelHmixDel4Init(err)
+
+ integer, intent(out) :: err
+
+ !--------------------------------------------------------------------
+ !
+ ! set some local module variables based on input config choices
+ !
+ !--------------------------------------------------------------------
+
+ err = 0
+
+ hmixDel4On = .false.
+
+ if ( config_h_mom_eddy_visc4 > 0.0 ) then
+ hmixDel4On = .true.
+ eddyVisc4 = config_h_mom_eddy_visc4
+ if (config_visc_vorticity_term) then
+ viscVortCoef = 1.0
+ else
+ viscVortCoef = 0.0
+ endif
+
+ endif
+
+ !--------------------------------------------------------------------
+
+ end subroutine OcnVelHmixDel4Init
+
+!***********************************************************************
+
+end module OcnVelHmixDel4
+
+!||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
Added: branches/ocean_projects/performance/src/core_ocean/module_OcnVelPressureGrad.F
===================================================================
--- branches/ocean_projects/performance/src/core_ocean/module_OcnVelPressureGrad.F (rev 0)
+++ branches/ocean_projects/performance/src/core_ocean/module_OcnVelPressureGrad.F 2011-09-15 21:59:53 UTC (rev 1002)
@@ -0,0 +1,196 @@
+!|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
+!
+! OcnVelPressureGrad
+!
+!> \brief MPAS ocean pressure gradient module
+!> \author Doug Jacobsen
+!> \date 15 September 2011
+!> \version SVN:$Id:$
+!> \details
+!> This module contains the routine for computing
+!> tendencie from the horizontal pressure gradient.
+!>
+!
+!-----------------------------------------------------------------------
+
+module OcnVelPressureGrad
+
+ use grid_types
+ use configure
+
+ implicit none
+ private
+ save
+
+ !--------------------------------------------------------------------
+ !
+ ! Public parameters
+ !
+ !--------------------------------------------------------------------
+
+ !--------------------------------------------------------------------
+ !
+ ! Public member functions
+ !
+ !--------------------------------------------------------------------
+
+ public :: OcnVelPressureGradTend, &
+ OcnVelPressureGradInit
+
+ !--------------------------------------------------------------------
+ !
+ ! Private module variables
+ !
+ !--------------------------------------------------------------------
+
+ real (kind=RKIND) :: rho0Inv
+
+
+!***********************************************************************
+
+contains
+
+!***********************************************************************
+!
+! routine OcnVelPressureGradTend
+!
+!> \brief Computes tendency term for horizontal pressure gradient
+!> \author Doug Jacobsen
+!> \date 15 September 2011
+!> \version SVN:$Id$
+!> \details
+!> This routine computes the pressure gradient tendency for momentum
+!> based on current state.
+!
+!-----------------------------------------------------------------------
+
+ subroutine OcnVelPressureGradTend(grid, pressure, tend, err)
+
+ !-----------------------------------------------------------------
+ !
+ ! input variables
+ !
+ !-----------------------------------------------------------------
+
+ real (kind=RKIND), dimension(:,:), intent(in) :: &
+ pressure !< Input: Pressure field or Mongomery potential
+
+ type (mesh_type), intent(in) :: &
+ grid !< Input: grid information
+
+ !-----------------------------------------------------------------
+ !
+ ! input/output variables
+ !
+ !-----------------------------------------------------------------
+
+ real (kind=RKIND), dimension(:,:), intent(inout) :: &
+ tend !< Input/Output: velocity tendency
+
+ !-----------------------------------------------------------------
+ !
+ ! output variables
+ !
+ !-----------------------------------------------------------------
+
+ integer, intent(out) :: err !< Output: error flag
+
+ !-----------------------------------------------------------------
+ !
+ ! local variables
+ !
+ !-----------------------------------------------------------------
+
+ integer :: nEdgesSolve, iEdge, k, cell1, cell2
+ integer, dimension(:), pointer :: maxLevelEdgeTop
+ integer, dimension(:,:), pointer :: cellsOnEdge
+
+ real (kind=RKIND), dimension(:), pointer :: dcEdge
+
+ err = 0
+
+ nEdgesSolve = grid % nEdgesSolve
+ maxLevelEdgeTop => grid % maxLevelEdgeTop % array
+ cellsOnEdge => grid % cellsOnEdge % array
+ dcEdge => grid % dcEdge % array
+
+ if (config_vert_grid_type.eq.'isopycnal') then
+ do iEdge=1,nEdgesSolve
+ cell1 = cellsOnEdge(1,iEdge)
+ cell2 = cellsOnEdge(2,iEdge)
+ do k=1,maxLevelEdgeTop(iEdge)
+ tend(k,iEdge) = tend(k,iEdge) &
+ - (pressure(k,cell2) - pressure(k,cell1))/dcEdge(iEdge)
+ end do
+ enddo
+ elseif (config_vert_grid_type.eq.'zlevel') then
+ do iEdge=1,nEdgesSolve
+ cell1 = cellsOnEdge(1,iEdge)
+ cell2 = cellsOnEdge(2,iEdge)
+ do k=1,maxLevelEdgeTop(iEdge)
+
+ tend(k,iEdge) = tend(k,iEdge) &
+ - rho0Inv*( pressure(k,cell2) &
+ - pressure(k,cell1) )/dcEdge(iEdge)
+ end do
+
+ enddo
+ endif
+
+
+ !--------------------------------------------------------------------
+
+ end subroutine OcnVelPressureGradTend
+
+!***********************************************************************
+!
+! routine OcnVelPressureGradInit
+!
+!> \brief Initializes ocean momentum horizontal pressure gradient
+!> \author Doug Jacobsen
+!> \date 15 September 2011
+!> \version SVN:$Id$
+!> \details
+!> This routine initializes parameters required for the computation of the
+!> horizontal pressure gradient.
+!
+!-----------------------------------------------------------------------
+
+
+ subroutine OcnVelPressureGradInit(err)
+
+ !--------------------------------------------------------------------
+
+
+ !-----------------------------------------------------------------
+ !
+ ! Output Variables
+ !
+ !-----------------------------------------------------------------
+
+ integer, intent(out) :: err !< Output: error flag
+
+
+ !-----------------------------------------------------------------
+ !
+ ! call individual init routines for each parameterization
+ !
+ !-----------------------------------------------------------------
+
+ err = 0
+
+ if (config_vert_grid_type.eq.'isopycnal') then
+ rho0Inv = 1.0
+ elseif (config_vert_grid_type.eq.'zlevel') then
+ rho0Inv = 1.0/config_rho0
+ end if
+
+ !--------------------------------------------------------------------
+
+ end subroutine OcnVelPressureGradInit
+
+!***********************************************************************
+
+end module OcnVelPressureGrad
+
+!|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
Added: branches/ocean_projects/performance/src/core_ocean/module_OcnVelVadv.F
===================================================================
--- branches/ocean_projects/performance/src/core_ocean/module_OcnVelVadv.F (rev 0)
+++ branches/ocean_projects/performance/src/core_ocean/module_OcnVelVadv.F 2011-09-15 21:59:53 UTC (rev 1002)
@@ -0,0 +1,193 @@
+!|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
+!
+! OcnVelVadv
+!
+!> \brief MPAS ocean vertical advection
+!> \author Doug Jacobsen
+!> \date 15 September 2011
+!> \version SVN:$Id:$
+!> \details
+!> This module contains the routine for computing
+!> tendencies for vertical advection.
+!>
+!
+!-----------------------------------------------------------------------
+
+module OcnVelVadv
+
+ use grid_types
+ use configure
+
+ implicit none
+ private
+ save
+
+ !--------------------------------------------------------------------
+ !
+ ! Public parameters
+ !
+ !--------------------------------------------------------------------
+
+ !--------------------------------------------------------------------
+ !
+ ! Public member functions
+ !
+ !--------------------------------------------------------------------
+
+ public :: OcnVelVadvTend, &
+ OcnVelVadvInit
+
+ !--------------------------------------------------------------------
+ !
+ ! Private module variables
+ !
+ !--------------------------------------------------------------------
+
+ logical :: velVadvOn
+
+
+!***********************************************************************
+
+contains
+
+!***********************************************************************
+!
+! routine OcnVelVadvTend
+!
+!> \brief Computes tendency term for vertical advection
+!> \author Doug Jacobsen
+!> \date 15 September 2011
+!> \version SVN:$Id$
+!> \details
+!> This routine computes the vertical advection tendency for momentum
+!> based on current state.
+!
+!-----------------------------------------------------------------------
+
+ subroutine OcnVelVadvTend(grid, u, wTop, tend, err)
+
+ !-----------------------------------------------------------------
+ !
+ ! input variables
+ !
+ !-----------------------------------------------------------------
+
+ real (kind=RKIND), dimension(:,:), intent(in) :: &
+ u, wTop
+
+ type (mesh_type), intent(in) :: &
+ grid !< Input: grid information
+
+ !-----------------------------------------------------------------
+ !
+ ! input/output variables
+ !
+ !-----------------------------------------------------------------
+
+ real (kind=RKIND), dimension(:,:), intent(inout) :: &
+ tend !< Input/Output: velocity tendency
+
+ !-----------------------------------------------------------------
+ !
+ ! output variables
+ !
+ !-----------------------------------------------------------------
+
+ integer, intent(out) :: err !< Output: error flag
+
+ !-----------------------------------------------------------------
+ !
+ ! local variables
+ !
+ !-----------------------------------------------------------------
+
+ integer :: iEdge, nEdgesSolve, cell1, cell2, k
+ integer :: nVertLevels
+ integer, dimension(:), pointer :: maxLevelEdgeTop
+ integer, dimension(:,:), pointer :: cellsOnEdge
+
+ real :: wTopEdge
+ real, dimension(:), allocatable :: w_dudzTopEdge
+ real, dimension(:), pointer :: zMidZLevel
+
+ if(.not.velVadvOn) return
+
+ err = 0
+
+ nVertLevels = grid % nVertLevels
+ nEdgesSolve = grid % nEdgesSolve
+ maxLevelEdgeTop => grid % maxLevelEdgeTop % array
+ cellsOnEdge => grid % cellsOnEdge % array
+ zMidZLevel => grid % zMidZLevel % array
+
+ allocate(w_dudzTopEdge(nVertLevels+1))
+ w_dudzTopEdge(1) = 0.0
+ do iEdge=1,nEdgesSolve
+ cell1 = cellsOnEdge(1,iEdge)
+ cell2 = cellsOnEdge(2,iEdge)
+
+ do k=2,maxLevelEdgeTop(iEdge)
+ ! Average w from cell center to edge
+ wTopEdge = 0.5*(wTop(k,cell1)+wTop(k,cell2))
+
+ ! compute dudz at vertical interface with first order derivative.
+ w_dudzTopEdge(k) = wTopEdge * (u(k-1,iEdge)-u(k,iEdge)) &
+ / (zMidZLevel(k-1) - zMidZLevel(k))
+ end do
+ w_dudzTopEdge(maxLevelEdgeTop(iEdge)+1) = 0.0
+ ! Average w*du/dz from vertical interface to vertical middle of cell
+ do k=1,maxLevelEdgeTop(iEdge)
+
+ tend(k,iEdge) = tend(k,iEdge) &
+ - 0.5 * (w_dudzTopEdge(k) + w_dudzTopEdge(k+1))
+ enddo
+ enddo
+ deallocate(w_dudzTopEdge)
+
+ !--------------------------------------------------------------------
+
+ end subroutine OcnVelVadvTend
+
+!***********************************************************************
+!
+! routine OcnVelVadvInit
+!
+!> \brief Initializes ocean momentum vertical advection
+!> \author Doug Jacobsen
+!> \date 15 September 2011
+!> \version SVN:$Id$
+!> \details
+!> This routine initializes a variety of quantities related to
+!> vertical velocity advection in the ocean.
+!
+!-----------------------------------------------------------------------
+
+
+ subroutine OcnVelVadvInit(err)
+
+ !--------------------------------------------------------------------
+
+ !-----------------------------------------------------------------
+ !
+ ! Output variables
+ !
+ !-----------------------------------------------------------------
+
+ integer, intent(out) :: err
+
+ err = 0
+ velVadvOn = .false.
+
+ if (config_vert_grid_type.eq.'zlevel') then
+ velVadvOn = .true.
+ end if
+
+ !--------------------------------------------------------------------
+
+ end subroutine OcnVelVadvInit
+
+!***********************************************************************
+
+end module OcnVelVadv
+
+!|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
Modified: branches/ocean_projects/performance/src/core_ocean/module_mpas_core.F
===================================================================
--- branches/ocean_projects/performance/src/core_ocean/module_mpas_core.F 2011-09-15 17:43:05 UTC (rev 1001)
+++ branches/ocean_projects/performance/src/core_ocean/module_mpas_core.F 2011-09-15 21:59:53 UTC (rev 1002)
@@ -5,6 +5,10 @@
use dmpar
use test_cases
+ use OcnVelPressureGrad
+ use OcnVelVadv
+ use OcnVelHmix
+
type (io_output_object) :: restart_obj
integer :: restart_frame
@@ -32,6 +36,8 @@
type (block_type), pointer :: block
type (dm_info) :: dminfo
+ integer :: err
+
if (.not. config_do_restart) call setup_sw_test_case(domain)
call compute_maxLevel(domain)
@@ -69,6 +75,10 @@
block => block % next
end do
+ call OcnVelPressureGradInit(err)
+ call OcnVelVadvInit(err)
+ call OcnVelHmixInit(err)
+
! mrp 100316 In order for this to work, we need to pass domain % dminfo as an
! input arguement into mpas_init. Ask about that later. For now, there will be
! no initial statistics write.
Modified: branches/ocean_projects/performance/src/core_ocean/module_time_integration.F
===================================================================
--- branches/ocean_projects/performance/src/core_ocean/module_time_integration.F 2011-09-15 17:43:05 UTC (rev 1001)
+++ branches/ocean_projects/performance/src/core_ocean/module_time_integration.F 2011-09-15 21:59:53 UTC (rev 1002)
@@ -8,6 +8,11 @@
use spline_interpolation
use timer
+ use OcnVelCoriolis
+ use OcnVelPressureGrad
+ use OcnVelVadv
+ use OcnVelHmix
+
contains
subroutine timestep(domain, dt, timeStamp)!{{{
@@ -1727,7 +1732,7 @@
integer :: iEdge, iCell, iVertex, k, cell1, cell2, &
vertex1, vertex2, eoe, i, j
- integer :: nCells, nEdges, nVertices, nVertLevels, nEdgesSolve
+ integer :: nCells, nEdges, nVertices, nVertLevels, nEdgesSolve, err
real (kind=RKIND) :: flux, vorticity_abs, h_vertex, workpv, q, &
upstream_bias, wTopEdge, rho0Inv, r, visc_vorticity_coef
real (kind=RKIND), dimension(:), pointer :: &
@@ -1820,85 +1825,31 @@
!
call timer_start("compute_tend_u-coriolis")
- do iEdge=1,grid % nEdgesSolve
- cell1 = cellsOnEdge(1,iEdge)
- cell2 = cellsOnEdge(2,iEdge)
- do k=1,maxLevelEdgeTop(iEdge)
+ call OcnVelCoriolisTend(grid, pv_edge, h_edge, u, ke, tend_u, err)
- q = 0.0
- do j = 1,nEdgesOnEdge(iEdge)
- eoe = edgesOnEdge(j,iEdge)
- workpv = 0.5 * (pv_edge(k,iEdge) + pv_edge(k,eoe))
- q = q + weightsOnEdge(j,iEdge) * u(k,eoe) * workpv * h_edge(k,eoe)
- end do
-
- tend_u(k,iEdge) = tend_u(k,iEdge) &
- + q &
- - ( ke(k,cell2) - ke(k,cell1) ) / dcEdge(iEdge)
-
- end do
- end do
call timer_stop("compute_tend_u-coriolis")
!
! velocity tendency: vertical advection term -w du/dz
!
call timer_start("compute_tend_u-vert adv")
- if (config_vert_grid_type.eq.'zlevel') then
- allocate(w_dudzTopEdge(nVertLevels+1))
- w_dudzTopEdge(1) = 0.0
- do iEdge=1,grid % nEdgesSolve
- cell1 = cellsOnEdge(1,iEdge)
- cell2 = cellsOnEdge(2,iEdge)
- do k=2,maxLevelEdgeTop(iEdge)
- ! Average w from cell center to edge
- wTopEdge = 0.5*(wTop(k,cell1)+wTop(k,cell2))
+ call OcnVelVadvTend(grid, u, wTop, tend_u, err)
- ! compute dudz at vertical interface with first order derivative.
- w_dudzTopEdge(k) = wTopEdge * (u(k-1,iEdge)-u(k,iEdge)) &
- / (zMidZLevel(k-1) - zMidZLevel(k))
- end do
- w_dudzTopEdge(maxLevelEdgeTop(iEdge)+1) = 0.0
- ! Average w*du/dz from vertical interface to vertical middle of cell
- do k=1,maxLevelEdgeTop(iEdge)
-
- tend_u(k,iEdge) = tend_u(k,iEdge) &
- - 0.5 * (w_dudzTopEdge(k) + w_dudzTopEdge(k+1))
- enddo
- enddo
- deallocate(w_dudzTopEdge)
- endif
call timer_stop("compute_tend_u-vert adv")
!
! velocity tendency: pressure gradient
!
call timer_start("compute_tend_u-pressure grad")
- rho0Inv = 1.0/config_rho0
+
if (config_vert_grid_type.eq.'isopycnal') then
- do iEdge=1,grid % nEdgesSolve
- cell1 = cellsOnEdge(1,iEdge)
- cell2 = cellsOnEdge(2,iEdge)
- do k=1,maxLevelEdgeTop(iEdge)
- tend_u(k,iEdge) = tend_u(k,iEdge) &
- - (MontPot(k,cell2) - MontPot(k,cell1))/dcEdge(iEdge)
- end do
- enddo
+ call OcnVelPressureGradTend(grid, MontPot, tend_u, err)
elseif (config_vert_grid_type.eq.'zlevel') then
- do iEdge=1,grid % nEdgesSolve
- cell1 = cellsOnEdge(1,iEdge)
- cell2 = cellsOnEdge(2,iEdge)
- do k=1,maxLevelEdgeTop(iEdge)
+ call OcnVelPressureGradTend(grid, pressure, tend_u, err)
+ end if
- tend_u(k,iEdge) = tend_u(k,iEdge) &
- - rho0Inv*( pressure(k,cell2) &
- - pressure(k,cell1) )/dcEdge(iEdge)
- end do
-
- enddo
- endif
call timer_stop("compute_tend_u-pressure grad")
!
@@ -1907,142 +1858,9 @@
! strictly only valid for config_h_mom_eddy_visc2 == constant
!
call timer_start("compute_tend_u-horiz mix")
- if (config_visc_vorticity_term) then
- visc_vorticity_coef = 1.0
- else
- visc_vorticity_coef = 0.0
- endif
-
- if ( config_h_mom_eddy_visc2 > 0.0 ) then
- call timer_start("compute_tend_u-horiz mix-del2")
- do iEdge=1,grid % nEdgesSolve
- cell1 = cellsOnEdge(1,iEdge)
- cell2 = cellsOnEdge(2,iEdge)
- vertex1 = verticesOnEdge(1,iEdge)
- vertex2 = verticesOnEdge(2,iEdge)
- do k=1,maxLevelEdgeTop(iEdge)
+ call OcnVelHmixTend(grid, divergence, vorticity, tend_u, err)
- ! Here -( vorticity(k,vertex2) - vorticity(k,vertex1) ) / dvEdge(iEdge)
- ! is - </font>
<font color="red">abla vorticity pointing from vertex 2 to vertex 1, or equivalently
- ! + k \times </font>
<font color="red">abla vorticity pointing from cell1 to cell2.
-
- u_diffusion = ( divergence(k,cell2) - divergence(k,cell1) ) / dcEdge(iEdge) &
- -visc_vorticity_coef &
- *( vorticity(k,vertex2) - vorticity(k,vertex1) ) / dvEdge(iEdge)
-
- u_diffusion = meshScalingDel2(iEdge) * config_h_mom_eddy_visc2 * u_diffusion
-
- tend_u(k,iEdge) = tend_u(k,iEdge) + u_diffusion
-
- end do
- end do
- call timer_stop("compute_tend_u-horiz mix-del2")
- end if
-
- !
- ! velocity tendency: del4 dissipation, -</font>
<font color="black">u_4 </font>
<font color="red">abla^4 u
- ! computed as </font>
<font color="black">abla^2 u = </font>
<font color="black">abla divergence + k \times </font>
<font color="red">abla vorticity
- ! applied recursively.
- ! strictly only valid for config_h_mom_eddy_visc4 == constant
- !
- if ( config_h_mom_eddy_visc4 > 0.0 ) then
- call timer_start("compute_tend_u-horiz mix-del4")
-
- allocate(delsq_divergence(nVertLevels, nCells+1))
- allocate(delsq_u(nVertLevels, nEdges+1))
- allocate(delsq_circulation(nVertLevels, nVertices+1))
- allocate(delsq_vorticity(nVertLevels, nVertices+1))
-
- delsq_u(:,:) = 0.0
-
- ! Compute </font>
<font color="black">abla^2 u = </font>
<font color="black">abla divergence + k \times </font>
<font color="red">abla vorticity
- do iEdge=1,grid % nEdges
- cell1 = cellsOnEdge(1,iEdge)
- cell2 = cellsOnEdge(2,iEdge)
- vertex1 = verticesOnEdge(1,iEdge)
- vertex2 = verticesOnEdge(2,iEdge)
-
- do k=1,maxLevelEdgeTop(iEdge)
-
- delsq_u(k,iEdge) = &
- ( divergence(k,cell2) - divergence(k,cell1) ) / dcEdge(iEdge) &
- -visc_vorticity_coef &
- *( vorticity(k,vertex2) - vorticity(k,vertex1)) / dvEdge(iEdge)
-
- end do
- end do
-
- ! vorticity using </font>
<font color="red">abla^2 u
- delsq_circulation(:,:) = 0.0
- do iEdge=1,nEdges
- vertex1 = verticesOnEdge(1,iEdge)
- vertex2 = verticesOnEdge(2,iEdge)
- do k=1,maxLevelEdgeTop(iEdge)
- delsq_circulation(k,vertex1) = delsq_circulation(k,vertex1) &
- - dcEdge(iEdge) * delsq_u(k,iEdge)
- delsq_circulation(k,vertex2) = delsq_circulation(k,vertex2) &
- + dcEdge(iEdge) * delsq_u(k,iEdge)
- end do
- end do
- do iVertex=1,nVertices
- r = 1.0 / areaTriangle(iVertex)
- do k=1,maxLevelVertexBot(iVertex)
- delsq_vorticity(k,iVertex) = delsq_circulation(k,iVertex) * r
- end do
- end do
-
- ! Divergence using </font>
<font color="red">abla^2 u
- delsq_divergence(:,:) = 0.0
- do iEdge=1,nEdges
- cell1 = cellsOnEdge(1,iEdge)
- cell2 = cellsOnEdge(2,iEdge)
- do k=1,maxLevelEdgeTop(iEdge)
- delsq_divergence(k,cell1) = delsq_divergence(k,cell1) &
- + delsq_u(k,iEdge)*dvEdge(iEdge)
- delsq_divergence(k,cell2) = delsq_divergence(k,cell2) &
- - delsq_u(k,iEdge)*dvEdge(iEdge)
- end do
- end do
- do iCell = 1,nCells
- r = 1.0 / areaCell(iCell)
- do k = 1,maxLevelCell(iCell)
- delsq_divergence(k,iCell) = delsq_divergence(k,iCell) * r
- end do
- end do
-
- ! Compute - \kappa </font>
<font color="red">abla^4 u
- ! as </font>
<font color="black">abla div(</font>
<font color="black">abla^2 u) + k \times </font>
<font color="black">abla ( k \cross curl(</font>
<font color="red">abla^2 u) )
- do iEdge=1,grid % nEdgesSolve
- cell1 = cellsOnEdge(1,iEdge)
- cell2 = cellsOnEdge(2,iEdge)
- vertex1 = verticesOnEdge(1,iEdge)
- vertex2 = verticesOnEdge(2,iEdge)
-
- do k=1,maxLevelEdgeTop(iEdge)
- delsq_u(k,iEdge) = &
- ( divergence(k,cell2) - divergence(k,cell1) ) / dcEdge(iEdge) &
- -( vorticity(k,vertex2) - vorticity(k,vertex1)) / dvEdge(iEdge)
-
- u_diffusion = ( delsq_divergence(k,cell2) &
- - delsq_divergence(k,cell1) ) / dcEdge(iEdge) &
- -visc_vorticity_coef &
- *( delsq_vorticity(k,vertex2) &
- - delsq_vorticity(k,vertex1) ) / dvEdge(iEdge)
-
- u_diffusion = meshScalingDel4(iEdge) * config_h_mom_eddy_visc4 * u_diffusion
-
- tend_u(k,iEdge) = tend_u(k,iEdge) - u_diffusion
- end do
- end do
-
- deallocate(delsq_divergence)
- deallocate(delsq_u)
- deallocate(delsq_circulation)
- deallocate(delsq_vorticity)
-
- call timer_stop("compute_tend_u-horiz mix-del4")
- end if
call timer_stop("compute_tend_u-horiz mix")
!
</font>
</pre>