<p><b>pwjones@lanl.gov</b> 2011-08-29 10:48:20 -0600 (Mon, 29 Aug 2011)</p><p><br>
changes to horizontal velocity mixing for the pointer hiding performance improvement and increased modularity for future additions<br>
</p><hr noshade><pre><font color="gray">Modified: branches/ocean_projects/performance/mpas/namelist.input.ocean
===================================================================
--- branches/ocean_projects/performance/mpas/namelist.input.ocean 2011-08-26 21:10:30 UTC (rev 957)
+++ branches/ocean_projects/performance/mpas/namelist.input.ocean 2011-08-29 16:48:20 UTC (rev 958)
@@ -24,8 +24,8 @@
config_rho0 = 1015
/
&hmix
- config_h_mom_eddy_visc2 = 0.0
- config_h_mom_eddy_visc4 = 5.0e8
+ config_OcnHmixVelDel2Visc = 0.0
+ config_OcnHmixVelDel4Visc = 5.0e8
config_OcnHmixTracerDel2Diff = 10.0
config_OcnHmixTracerDel4Diff = 0.0
/
Modified: branches/ocean_projects/performance/mpas/namelist.input.sw
===================================================================
--- branches/ocean_projects/performance/mpas/namelist.input.sw 2011-08-26 21:10:30 UTC (rev 957)
+++ branches/ocean_projects/performance/mpas/namelist.input.sw 2011-08-29 16:48:20 UTC (rev 958)
@@ -24,8 +24,8 @@
config_rho0 = 1000
/
&hmix
- config_h_mom_eddy_visc2 = 1.0e5
- config_h_mom_eddy_visc4 = 0.0
+ config_OcnHmixVelDel2Visc = 1.0e5
+ config_OcnHmixVelDel4Visc = 0.0
config_OcnHmixTracerDel2Diff = 1.0e4
config_OcnHmixTracerDel4Diff = 0.0
/
Modified: branches/ocean_projects/performance/mpas/src/core_ocean/Makefile
===================================================================
--- branches/ocean_projects/performance/mpas/src/core_ocean/Makefile 2011-08-26 21:10:30 UTC (rev 957)
+++ branches/ocean_projects/performance/mpas/src/core_ocean/Makefile 2011-08-29 16:48:20 UTC (rev 958)
@@ -6,10 +6,13 @@
module_OcnHmixTracer.o \
module_OcnHmixTracerDel2.o \
module_OcnHmixTracerDel4.o \
- module_OcnVmixTracer.o \
- module_OcnVmixTracerConst.o \
- module_OcnVmixTracerTanh.o \
- module_OcnCoriolis.o \
+ module_OcnHmixVel.o \
+ module_OcnHmixVelDel2.o \
+ module_OcnHmixVelDel4.o \
+ module_OcnVmixTracer.o \
+ module_OcnVmixTracerConst.o \
+ module_OcnVmixTracerTanh.o \
+ module_OcnCoriolis.o \
module_time_integration.o \
module_global_diagnostics.o
@@ -36,11 +39,17 @@
module_OcnHmixTracer.o: module_OcnHmixTracerDel2.o module_OcnHmixTracerDel4.o
+module_OcnHmixVelDel2.o:
+
+module_OcnHmixVelDel4.o:
+
+module_OcnHmixVel.o: module_OcnHmixVelDel2.o module_OcnHmixVelDel4.o
+
module_time_integration.o:
module_global_diagnostics.o:
-module_mpas_core.o: module_advection.o module_global_diagnostics.o module_test_cases.o module_OcnCoriolis.o module_OcnVmixTracer.o module_OcnVmixTracerConst.o module_OcnVmixTracerTanh.o module_OcnHmixTracer.o module_OcnHmixTracerDel2.o module_OcnHmixTracerDel4.o module_time_integration.o
+module_mpas_core.o: module_advection.o module_global_diagnostics.o module_test_cases.o module_OcnCoriolis.o module_OcnVmixTracer.o module_OcnVmixTracerConst.o module_OcnVmixTracerTanh.o module_OcnHmixTracer.o module_OcnHmixTracerDel2.o module_OcnHmixTracerDel4.o module_OcnHmixVel.o module_OcnHmixVelDel2.o module_OcnHmixVelDel4.o module_time_integration.o
clean:
$(RM) *.o *.mod *.f90 libdycore.a
Modified: branches/ocean_projects/performance/mpas/src/core_ocean/Registry
===================================================================
--- branches/ocean_projects/performance/mpas/src/core_ocean/Registry 2011-08-26 21:10:30 UTC (rev 957)
+++ branches/ocean_projects/performance/mpas/src/core_ocean/Registry 2011-08-29 16:48:20 UTC (rev 958)
@@ -16,8 +16,8 @@
namelist real restart config_restart_time 172800.0
namelist character grid config_vert_grid_type isopycnal
namelist real grid config_rho0 1028
-namelist real hmix config_h_mom_eddy_visc2 0.0
-namelist real hmix config_h_mom_eddy_visc4 0.0
+namelist real hmix config_OcnHmixVelDel2Visc 0.0
+namelist real hmix config_OcnHmixVelDel4Visc 0.0
namelist real hmix config_OcnHmixTracerDel2Diff 0.0
namelist real hmix config_OcnHmixTracerDel4Diff 0.0
namelist real hmix config_apvm_upwinding 0.5
Added: branches/ocean_projects/performance/mpas/src/core_ocean/module_OcnHmixVel.F
===================================================================
--- branches/ocean_projects/performance/mpas/src/core_ocean/module_OcnHmixVel.F (rev 0)
+++ branches/ocean_projects/performance/mpas/src/core_ocean/module_OcnHmixVel.F 2011-08-29 16:48:20 UTC (rev 958)
@@ -0,0 +1,164 @@
+!|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
+!
+! OcnHmixVel
+!
+!> \brief MPAS ocean horizontal momentum mixing driver
+!> \author Phil Jones
+!> \date 26 July 2011
+!> \version SVN:$Id:$
+!> \details
+!> This module contains the main driver routine for computing
+!> horizontal mixing tendencies. It primarily decides which mixing
+!> parameterizations are being used and calls them individually.
+!> Detailed mixing parameterizations are contained in their own
+!> modules.
+!
+!-----------------------------------------------------------------------
+
+module OcnHmixVel
+
+ use grid_types
+ use configure
+ use OcnHmixVelDel2
+ use OcnHmixVelDel4
+
+ implicit none
+ private
+ save
+
+ !--------------------------------------------------------------------
+ !
+ ! Public parameters
+ !
+ !--------------------------------------------------------------------
+
+ !--------------------------------------------------------------------
+ !
+ ! Public member functions
+ !
+ !--------------------------------------------------------------------
+
+ public :: OcnHmixVelTend, &
+ OcnHmixVelInit
+
+ !--------------------------------------------------------------------
+ !
+ ! Private module variables
+ !
+ !--------------------------------------------------------------------
+
+
+!***********************************************************************
+
+contains
+
+!***********************************************************************
+!
+! routine OcnHmixVelTend
+!
+!> \brief Computes tendency term for horizontal momentum mixing
+!> \author Phil Jones
+!> \date 26 July 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 OcnHmixVelTend(tend, divergence, vorticity, grid)
+
+ !-----------------------------------------------------------------
+ !
+ ! 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
+ !
+ !-----------------------------------------------------------------
+
+ !-----------------------------------------------------------------
+ !
+ ! local variables
+ !
+ !-----------------------------------------------------------------
+
+ !-----------------------------------------------------------------
+ !
+ ! call relevant routines for computing tendencies
+ ! note that the user can choose multiple options and the
+ ! tendencies will be added together
+ !
+ !-----------------------------------------------------------------
+
+ call OcnHmixVelDel2Tend(tend, divergence, vorticity, grid)
+ call OcnHmixVelDel4Tend(tend, divergence, vorticity, grid)
+
+ !--------------------------------------------------------------------
+
+ end subroutine OcnHmixVelTend
+
+!***********************************************************************
+!
+! routine OcnHmixVelInit
+!
+!> \brief Initializes ocean momentum horizontal mixing quantities
+!> \author Phil Jones
+!> \date 26 July 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 OcnHmixVelInit
+
+ !--------------------------------------------------------------------
+
+ !-----------------------------------------------------------------
+ !
+ ! call individual init routines for each parameterization
+ !
+ !-----------------------------------------------------------------
+
+ call OcnHmixVelDel2Init
+ call OcnHmixVelDel4Init
+
+ !--------------------------------------------------------------------
+
+ end subroutine OcnHmixVelInit
+
+!***********************************************************************
+
+end module OcnHmixVel
+
+!|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
Added: branches/ocean_projects/performance/mpas/src/core_ocean/module_OcnHmixVelDel2.F
===================================================================
--- branches/ocean_projects/performance/mpas/src/core_ocean/module_OcnHmixVelDel2.F (rev 0)
+++ branches/ocean_projects/performance/mpas/src/core_ocean/module_OcnHmixVelDel2.F 2011-08-29 16:48:20 UTC (rev 958)
@@ -0,0 +1,215 @@
+!|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
+!
+! OcnHMixVelDel2
+!
+!> \brief Ocean horizontal mixing - Laplacian parameterization
+!> \author Phil Jones
+!> \date 22 August 2011
+!> \version SVN:$Id:$
+!> \details
+!> This module contains routines for computing horizontal mixing
+!> tendencies using a Laplacian formulation.
+!
+!-----------------------------------------------------------------------
+
+module OcnHMixVelDel2
+
+ use grid_types
+ use configure
+
+ implicit none
+ private
+ save
+
+ !--------------------------------------------------------------------
+ !
+ ! Public parameters
+ !
+ !--------------------------------------------------------------------
+
+ !--------------------------------------------------------------------
+ !
+ ! Public member functions
+ !
+ !--------------------------------------------------------------------
+
+ public :: OcnHMixVelDel2Tend, &
+ OcnHMixVelDel2Init
+
+ !-------------------------------------------------------------------
+ !
+ ! Private module variables
+ !
+ !--------------------------------------------------------------------
+
+ logical :: &
+ hmixDel2On !< local flag to determine whether del2 chosen
+
+ real (kind=RKIND) :: &
+ eddyVisc2 !< base eddy diffusivity for Laplacian
+
+!***********************************************************************
+
+contains
+
+!***********************************************************************
+!
+! routine OcnHMixVelDel2Tend
+!
+!> \brief Computes tendency term for Laplacian horizontal momentum mixing
+!> \author Phil Jones
+!> \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 OcnHMixVelDel2Tend(tend, divergence, vorticity, grid)
+
+ !-----------------------------------------------------------------
+ !
+ ! 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
+
+ !-----------------------------------------------------------------
+ !
+ ! local variables
+ !
+ !-----------------------------------------------------------------
+
+ integer :: &
+ iEdge, k, &! loop counters
+ cell1, cell2, &! cell indices
+ vertex1, vertex2, &! vertex indices
+ nEdgesSolve ! array sizes
+
+ real (kind=RKIND), dimension(:), pointer :: &
+ dvEdge, dcEdge
+
+ integer, dimension(:), pointer :: maxLevelEdgeTop
+
+ integer, dimension(:,:), pointer :: &
+ cellsOnEdge, &! cell index on edges
+ verticesOnEdge ! vertex index on edges
+
+ !-----------------------------------------------------------------
+ !
+ ! exit if this mixing is not selected
+ !
+ !-----------------------------------------------------------------
+
+ if (.not. hmixDel2On) return
+
+ !-----------------------------------------------------------------
+ !
+ ! initialize some arrays and sizes
+ !
+ !-----------------------------------------------------------------
+
+ cellsOnEdge => grid % cellsOnEdge % array
+ verticesOnEdge => grid % verticesOnEdge % array
+ dvEdge => grid % dvEdge % array
+ dcEdge => grid % dcEdge % array
+ maxLevelEdgeTop => grid % maxLevelEdgeTop % array
+
+ nEdgesSolve = grid % nEdgesSolve
+
+ !-----------------------------------------------------------------
+ !
+ ! del2 horizontal momentum diffusion
+ !
+ ! Here -( vort(k,vertex2) - vort(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.
+ !
+ !-----------------------------------------------------------------
+
+ 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)
+
+ tend(k,iEdge) = tend(k,iEdge) + eddyVisc2 *( &
+ ( divergence(k,cell2 ) - &
+ divergence(k,cell1 ) ) / dcEdge(iEdge) &
+ -( vorticity (k,vertex2) - &
+ vorticity (k,vertex1) ) / dvEdge(iEdge))
+
+ end do
+ end do
+
+ !--------------------------------------------------------------------
+
+ end subroutine OcnHMixVelDel2Tend
+
+!***********************************************************************
+!
+! routine OcnHMixVelDel2Init
+!
+!> \brief Initializes ocean momentum Laplacian horizontal mixing
+!> \author Phil Jones
+!> \date 26 July 2011
+!> \version SVN:$Id$
+!> \details
+!> This routine initializes a variety of quantities related to
+!> Laplacian horizontal momentum mixing in the ocean.
+!
+!-----------------------------------------------------------------------
+
+ subroutine OcnHMixVelDel2Init
+
+ !--------------------------------------------------------------------
+ !
+ ! set some local module variables based on input config choices
+ !
+ !--------------------------------------------------------------------
+
+ hmixDel2On = .false.
+
+ if (config_OcnHmixVelDel2Visc > 0.0 ) then
+
+ hmixDel2On = .true.
+
+ eddyVisc2 = config_OcnHmixVelDel2Visc
+
+ endif
+
+ !--------------------------------------------------------------------
+
+ end subroutine OcnHMixVelDel2Init
+
+!***********************************************************************
+
+end module OcnHMixVelDel2
+
+!|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
Added: branches/ocean_projects/performance/mpas/src/core_ocean/module_OcnHmixVelDel4.F
===================================================================
--- branches/ocean_projects/performance/mpas/src/core_ocean/module_OcnHmixVelDel4.F (rev 0)
+++ branches/ocean_projects/performance/mpas/src/core_ocean/module_OcnHmixVelDel4.F 2011-08-29 16:48:20 UTC (rev 958)
@@ -0,0 +1,317 @@
+!|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
+!
+! OcnHMixVelDel4
+!
+!> \brief Ocean horizontal mixing - biharmonic parameterization
+!> \author Phil Jones
+!> \date 1 August 2011
+!> \version SVN:$Id:$
+!> \details
+!> This module contains routines and variables for computing
+!> horizontal mixing tendencies using a biharmonic formulation.
+!
+!-----------------------------------------------------------------------
+
+module OcnHMixVelDel4
+
+ use grid_types
+ use configure
+
+ implicit none
+ private
+ save
+
+ !--------------------------------------------------------------------
+ !
+ ! Public parameters
+ !
+ !--------------------------------------------------------------------
+
+ !--------------------------------------------------------------------
+ !
+ ! Public member functions
+ !
+ !--------------------------------------------------------------------
+
+ public :: OcnHMixVelDel4Tend, &
+ OcnHMixVelDel4Init
+
+ !--------------------------------------------------------------------
+ !
+ ! Private module variables
+ !
+ !--------------------------------------------------------------------
+
+ logical :: &
+ hmixDel4On !< local flag to determine whether del4 chosen
+
+ real (kind=RKIND) :: &
+ eddyVisc4 !< base eddy diffusivity for biharmonic
+
+
+!***********************************************************************
+
+contains
+
+!***********************************************************************
+!
+! routine OcnHMixVelDel4Tend
+!
+!> \brief Computes tendency term for biharmonic horizontal momentum mixing
+!> \author Phil Jones
+!> \date 22 August 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 OcnHMixVelDel4Tend(tend, divergence, vorticity, grid)
+
+ !-----------------------------------------------------------------
+ !
+ ! 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
+
+ !-----------------------------------------------------------------
+ !
+ ! local variables
+ !
+ !-----------------------------------------------------------------
+
+ integer :: &
+ i, k, iCell, iEdge, &! loop counters
+ cell1, cell2, &! cell indices
+ vertex1, vertex2, &! vertex indices
+ nEdges, nCells, nVertices, &! array sizes
+ nVertLevels ! array sizes
+
+ real (kind=RKIND) :: &
+ invArea1, invArea2, &! reciprocal of cell areas
+ flux ! flux across edge
+
+ real (kind=RKIND), dimension(:), pointer :: &
+ dvEdge, dcEdge, areaCell, &! grid quantities
+ areaTriangle
+
+ integer, dimension(:), pointer :: &
+ maxLevelCell, &! max vert level index for cells
+ maxLevelEdgeTop, &! max vert level index for top edge
+ maxLevelVertexBot ! max vert level index for bot vertex
+
+ integer, dimension(:,:), pointer :: &
+ cellsOnEdge, &! cell index on edges
+ verticesOnEdge ! vertex index on edges
+
+ real (kind=RKIND), dimension(:,:), allocatable :: &
+ delsq_u, &! temp for first del2 application
+ delsq_divergence, &! temp for first del2 application
+ delsq_vorticity ! temp for first del2 application
+
+ !-----------------------------------------------------------------
+ !
+ ! exit if this mixing is not selected
+ !
+ !-----------------------------------------------------------------
+
+ if (.not. hmixDel4On) return
+
+ !-----------------------------------------------------------------
+ !
+ ! initialize some arrays and sizes
+ !
+ !-----------------------------------------------------------------
+
+ areaCell => grid % areaCell % array
+ areaTriangle => grid % areaTriangle % array
+ cellsOnEdge => grid % cellsOnEdge % array
+ verticesOnEdge => grid % verticesOnEdge % array
+ dvEdge => grid % dvEdge % array
+ dcEdge => grid % dcEdge % array
+ maxLevelCell => grid % maxLevelCell % array
+ maxLevelEdgeTop => grid % maxLevelEdgeTop % array
+ maxLevelVertexBot => grid % maxLevelVertexBot % array
+
+ nEdges = grid % nEdges
+ nCells = grid % nCells
+ nVertices = grid % nVertices
+ nVertLevels = size(vorticity,dim=1)
+
+ allocate(delsq_u (nVertLevels, nEdges+1))
+ allocate(delsq_divergence (nVertLevels, nCells+1))
+ allocate(delsq_vorticity (nVertLevels, nVertices+1))
+
+ !-----------------------------------------------------------------
+ !
+ ! Compute </font>
<font color="black">abla^2 u = </font>
<font color="blue">abla divergence +
+ ! k \times </font>
<font color="blue">abla vorticity
+ !
+ !-----------------------------------------------------------------
+
+ do iEdge=1,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) &
+ -( vorticity(k,vertex2) - &
+ vorticity(k,vertex1)) / dvEdge(iEdge)
+
+ end do
+ do k=maxLevelEdgeTop(iEdge)+1,nVertLevels
+ delsq_u(:,:) = 0.d0
+ end do
+
+ end do
+
+ !***
+ !*** vorticity using </font>
<font color="blue">abla^2 u
+ !***
+
+ delsq_vorticity(:,:) = 0.0
+
+ do iEdge=1,nEdges
+
+ vertex1 = verticesOnEdge(1,iEdge)
+ vertex2 = verticesOnEdge(2,iEdge)
+ invArea1 = 1.d0 / areaTriangle(vertex1)
+ invArea2 = 1.d0 / areaTriangle(vertex2)
+
+ do k=1,maxLevelEdgeTop(iEdge)
+ delsq_vorticity(k,vertex1) = delsq_vorticity(k,vertex1) &
+ - dcEdge(iEdge) * delsq_u(k,iEdge) * invArea1
+ delsq_vorticity(k,vertex2) = delsq_vorticity(k,vertex2) &
+ + dcEdge(iEdge) * delsq_u(k,iEdge) * invArea2
+ end do
+ end do
+
+ !***
+ !*** Divergence using </font>
<font color="blue">abla^2 u
+ !***
+
+ delsq_divergence(:,:) = 0.d0
+
+ do iEdge=1,nEdges
+
+ cell1 = cellsOnEdge(1,iEdge)
+ cell2 = cellsOnEdge(2,iEdge)
+ invArea1 = 1.d0 / areaCell(cell1)
+ invArea2 = 1.d0 / areaCell(cell2)
+
+ do k=1,maxLevelEdgeTop(iEdge)
+ delsq_divergence(k,cell1) = delsq_divergence(k,cell1) &
+ + delsq_u(k,iEdge)*dvEdge(iEdge)*invArea1
+ delsq_divergence(k,cell2) = delsq_divergence(k,cell2) &
+ - delsq_u(k,iEdge)*dvEdge(iEdge)*invArea2
+ end do
+ end do
+
+ !-----------------------------------------------------------------
+ !
+ ! Compute \kappa </font>
<font color="blue">abla^4 u
+ ! as </font>
<font color="black">abla div(</font>
<font color="blue">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)
+
+ tend(k,iEdge) = tend(k,iEdge) - eddyVisc4 * ( &
+ ( delsq_divergence(k,cell2 ) &
+ - delsq_divergence(k,cell1 ) ) / dcEdge(iEdge) &
+ -( delsq_vorticity (k,vertex2) &
+ - delsq_vorticity (k,vertex1) ) / dvEdge(iEdge))
+
+ end do
+ end do
+
+ !***
+ !*** release storage
+ !***
+
+ deallocate(delsq_divergence)
+ deallocate(delsq_u)
+ deallocate(delsq_vorticity)
+
+ !--------------------------------------------------------------------
+
+ end subroutine OcnHMixVelDel4Tend
+
+!***********************************************************************
+!
+! routine OcnHMixVelDel4Init
+!
+!> \brief Initializes ocean momentum biharmonic horizontal mixing
+!> \author Phil Jones
+!> \date 22 August 2011
+!> \version SVN:$Id$
+!> \details
+!> This routine initializes a variety of quantities related to
+!> biharmonic horizontal tracer mixing in the ocean.
+!
+!-----------------------------------------------------------------------
+
+
+ subroutine OcnHMixVelDel4Init
+
+ !--------------------------------------------------------------------
+ !
+ ! set some local module variables based on input config choices
+ !
+ !--------------------------------------------------------------------
+
+ hmixDel4On = .false.
+
+ if (config_OcnHmixVelDel4Visc /= 0.0 ) then
+
+ hmixDel4On = .true.
+
+ eddyVisc4 = config_OcnHmixVelDel4Visc
+
+ endif
+
+ !--------------------------------------------------------------------
+
+ end subroutine OcnHMixVelDel4Init
+
+!***********************************************************************
+
+end module OcnHMixVelDel4
+
+!||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
Modified: branches/ocean_projects/performance/mpas/src/core_ocean/module_mpas_core.F
===================================================================
--- branches/ocean_projects/performance/mpas/src/core_ocean/module_mpas_core.F 2011-08-26 21:10:30 UTC (rev 957)
+++ branches/ocean_projects/performance/mpas/src/core_ocean/module_mpas_core.F 2011-08-29 16:48:20 UTC (rev 958)
@@ -3,6 +3,7 @@
use mpas_framework
use dmpar
use test_cases
+ use OcnHmixVel
use OcnHmixTracer
use OcnVmixTracer
@@ -57,6 +58,7 @@
!***
call OcnHmixTracerInit
+ call OcnHmixVelInit
call OcnVmixTracerInit(err)
if(err .ne. 0) then
Modified: branches/ocean_projects/performance/mpas/src/core_ocean/module_time_integration.F
===================================================================
--- branches/ocean_projects/performance/mpas/src/core_ocean/module_time_integration.F 2011-08-26 21:10:30 UTC (rev 957)
+++ branches/ocean_projects/performance/mpas/src/core_ocean/module_time_integration.F 2011-08-29 16:48:20 UTC (rev 958)
@@ -9,6 +9,7 @@
use vector_reconstruction
use spline_interpolation
use timer
+ use OcnHmixVel
use OcnHmixTracer
use OcnVmixTracer
use OcnCoriolis
@@ -133,7 +134,7 @@
block % mesh % nVertLevels, block % mesh % nEdges, &
block % parinfo % edgesToSend, block % parinfo % edgesToRecv)
- if (config_h_mom_eddy_visc4 > 0.0) then
+ if (config_OcnHmixVelDel4Visc > 0.0) then
call dmpar_exch_halo_field2dReal(domain % dminfo, block % state % time_levs(2) % state % divergence % array(:,:), &
block % mesh % nVertLevels, block % mesh % nCells, &
block % parinfo % cellsToSend, block % parinfo % cellsToRecv)
@@ -459,121 +460,12 @@
endif
call timer_stop("press grad")
+
call timer_start("vel dissipation")
- !
- ! velocity tendency: del2 dissipation, </font>
<font color="black">u_2 </font>
<font color="red">abla^2 u
- ! computed as </font>
<font color="black">u( </font>
<font color="black">abla divergence + k \times </font>
<font color="red">abla vorticity )
- ! strictly only valid for config_h_mom_eddy_visc2 == constant
- !
- if ( config_h_mom_eddy_visc2 > 0.0 ) then
- do iEdge=1,grid % nEdgesSolve
- cell1 = cellsOnEdge(1,iEdge)
- cell2 = cellsOnEdge(2,iEdge)
- vertex1 = verticesOnEdge(1,iEdge)
- vertex2 = verticesOnEdge(2,iEdge)
+ call OcnHmixVelTend(tend_u, divergence, vorticity, grid)
+ call timer_stop("vel dissipation")
- do k=1,maxLevelEdgeTop(iEdge)
-
- ! 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) &
- -( vorticity(k,vertex2) - vorticity(k,vertex1) ) / dvEdge(iEdge)
- u_diffusion = config_h_mom_eddy_visc2 * u_diffusion
-
- tend_u(k,iEdge) = tend_u(k,iEdge) + u_diffusion
-
- end do
- end do
- 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
-
- 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
- delsq_circulation(:,:) = 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) &
- -( vorticity(k,vertex2) - vorticity(k,vertex1)) / dvEdge(iEdge)
-
- !---- DWJ Testing Comments for loop merging
- end do
- end do
-
- ! vorticity using </font>
<font color="red">abla^2 u
- 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) / areaTriangle(vertex1)
- delsq_circulation(k,vertex2) = delsq_circulation(k,vertex2) &
- + dcEdge(iEdge) * delsq_u(k,iEdge) / areaTriangle(vertex2)
- 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) / areaCell(cell1)
- delsq_divergence(k,cell2) = delsq_divergence(k,cell2) &
- - delsq_u(k,iEdge)*dvEdge(iEdge) / areaCell(cell2)
- 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)
-
- u_diffusion = ( delsq_divergence(k,cell2) &
- - delsq_divergence(k,cell1) ) / dcEdge(iEdge) &
- -( delsq_vorticity(k,vertex2) &
- - delsq_vorticity(k,vertex1) ) / dvEdge(iEdge)
-
- tend_u(k,iEdge) = tend_u(k,iEdge) - config_h_mom_eddy_visc4 * u_diffusion
- end do
- end do
-
- deallocate(delsq_divergence)
- deallocate(delsq_u)
- deallocate(delsq_circulation)
- deallocate(delsq_vorticity)
-
- end if
-
- call timer_stop("vel dissipation")
- !
! velocity tendency: nonlinear Coriolis term and grad of kinetic energy
!
call timer_start("coriolis")
</font>
</pre>