<p><b>mpetersen@lanl.gov</b> 2012-10-22 13:36:32 -0600 (Mon, 22 Oct 2012)</p><p>trunk commit: Merge leith_mrp branch to trunk. Reviewed by Doug. I added a new module for the Leith closure. For non-Leith runs, I get bit-for-bit match with trunk. With Leith, get b-f-b match between 16 and 64 procs on the 120km global. Also tested with overflow domain. The Leith closure was originally written by Todd.<br>
</p><hr noshade><pre><font color="gray">Index: trunk/mpas
===================================================================
--- trunk/mpas 2012-10-22 19:14:22 UTC (rev 2241)
+++ trunk/mpas 2012-10-22 19:36:32 UTC (rev 2242)
Property changes on: trunk/mpas
___________________________________________________________________
Modified: svn:mergeinfo
## -6,6 +6,7 ##
/branches/ocean_projects/gmvar:1214-1514,1517-1738
/branches/ocean_projects/imp_vert_mix_error:1847-1887
/branches/ocean_projects/imp_vert_mix_mrp:754-986
+/branches/ocean_projects/leith_mrp:2182-2241
/branches/ocean_projects/monotonic_advection:1499-1640
/branches/ocean_projects/monthly_forcing:1810-1867
/branches/ocean_projects/option3_b4b_test:2201-2231
\ No newline at end of property
Modified: trunk/mpas/namelist.input.ocean
===================================================================
--- trunk/mpas/namelist.input.ocean 2012-10-22 19:14:22 UTC (rev 2241)
+++ trunk/mpas/namelist.input.ocean 2012-10-22 19:36:32 UTC (rev 2242)
@@ -62,6 +62,12 @@
config_h_tracer_eddy_diff2 = 1.0e5
config_h_tracer_eddy_diff4 = 0.0
/
+&hmix_leith
+ config_use_leith_del2 = .false.
+ config_leith_parameter = 1.0
+ config_leith_dx = 15000.0
+ config_leith_visc2_max = 2.5e3
+/
&vmix
config_vert_visc_type = 'const'
config_vert_diff_type = 'const'
Index: trunk/mpas/src/core_ocean
===================================================================
--- trunk/mpas/src/core_ocean 2012-10-22 19:14:22 UTC (rev 2241)
+++ trunk/mpas/src/core_ocean 2012-10-22 19:36:32 UTC (rev 2242)
Property changes on: trunk/mpas/src/core_ocean
___________________________________________________________________
Modified: svn:mergeinfo
## -6,6 +6,7 ##
/branches/ocean_projects/gmvar/src/core_ocean:1214-1514,1517-1738
/branches/ocean_projects/imp_vert_mix_error/src/core_ocean:1847-1887
/branches/ocean_projects/imp_vert_mix_mrp/src/core_ocean:754-986
+/branches/ocean_projects/leith_mrp/src/core_ocean:2182-2241
/branches/ocean_projects/monotonic_advection/src/core_ocean:1499-1640
/branches/ocean_projects/monthly_forcing/src/core_ocean:1810-1867
/branches/ocean_projects/option3_b4b_test/src/core_ocean:2201-2231
\ No newline at end of property
Modified: trunk/mpas/src/core_ocean/Makefile
===================================================================
--- trunk/mpas/src/core_ocean/Makefile 2012-10-22 19:14:22 UTC (rev 2241)
+++ trunk/mpas/src/core_ocean/Makefile 2012-10-22 19:36:32 UTC (rev 2242)
@@ -10,6 +10,7 @@
mpas_ocn_vel_vadv.o \
mpas_ocn_vel_hmix.o \
mpas_ocn_vel_hmix_del2.o \
+ mpas_ocn_vel_hmix_leith.o \
mpas_ocn_vel_hmix_del4.o \
mpas_ocn_vel_forcing.o \
mpas_ocn_vel_forcing_windstress.o \
@@ -87,10 +88,12 @@
mpas_ocn_vel_vadv.o:
-mpas_ocn_vel_hmix.o: mpas_ocn_vel_hmix_del2.o mpas_ocn_vel_hmix_del4.o
+mpas_ocn_vel_hmix.o: mpas_ocn_vel_hmix_del2.o mpas_ocn_vel_hmix_leith.o mpas_ocn_vel_hmix_del4.o
mpas_ocn_vel_hmix_del2.o:
+mpas_ocn_vel_hmix_leith.o:
+
mpas_ocn_vel_hmix_del4.o:
mpas_ocn_vel_forcing.o: mpas_ocn_vel_forcing_windstress.o mpas_ocn_vel_forcing_bottomdrag.o mpas_ocn_vel_forcing_rayleigh.o
@@ -179,6 +182,7 @@
mpas_ocn_vel_vadv.o \
mpas_ocn_vel_hmix.o \
mpas_ocn_vel_hmix_del2.o \
+ mpas_ocn_vel_hmix_leith.o \
mpas_ocn_vel_hmix_del4.o \
mpas_ocn_vel_forcing.o \
mpas_ocn_vel_forcing_windstress.o \
Modified: trunk/mpas/src/core_ocean/Registry
===================================================================
--- trunk/mpas/src/core_ocean/Registry 2012-10-22 19:14:22 UTC (rev 2241)
+++ trunk/mpas/src/core_ocean/Registry 2012-10-22 19:36:32 UTC (rev 2242)
@@ -63,6 +63,10 @@
namelist logical hmix config_rayleigh_friction false
namelist real hmix config_rayleigh_damping_coeff 0.0
namelist real hmix config_apvm_scale_factor 0.0
+namelist logical hmix_leith config_use_leith_del2 false
+namelist real hmix_leith config_leith_parameter 0.0
+namelist real hmix_leith config_leith_dx 0.0
+namelist real hmix_leith config_leith_visc2_max 1000000.0
namelist character vmix config_vert_visc_type const
namelist character vmix config_vert_diff_type const
namelist logical vmix config_implicit_vertical_mix true
@@ -281,6 +285,7 @@
var persistent real pressure ( nVertLevels nCells Time ) 2 - pressure state - -
var persistent real wTop ( nVertLevelsP1 nCells Time ) 2 - wTop state - -
var persistent real rhoDisplaced ( nVertLevels nCells Time ) 2 - rhoDisplaced state - -
+var persistent real viscosity ( nVertLevels nEdges Time ) 2 o viscosity state - -
% Other diagnostic variables: neither read nor written to any files
var persistent real vh ( nVertLevels nEdges Time ) 2 - vh state - -
Modified: trunk/mpas/src/core_ocean/mpas_ocn_tendency.F
===================================================================
--- trunk/mpas/src/core_ocean/mpas_ocn_tendency.F 2012-10-22 19:14:22 UTC (rev 2241)
+++ trunk/mpas/src/core_ocean/mpas_ocn_tendency.F 2012-10-22 19:36:32 UTC (rev 2242)
@@ -172,7 +172,7 @@
real (kind=RKIND), dimension(:,:), pointer :: &
h_edge, h, u, rho, zMid, pressure, &
- tend_u, circulation, vorticity, ke, ke_edge, Vor_edge, &
+ tend_u, circulation, vorticity, viscosity, ke, ke_edge, Vor_edge, &
MontPot, wTop, divergence, vertViscTopOfEdge
real (kind=RKIND), dimension(:,:), pointer :: u_src
@@ -186,6 +186,7 @@
wTop => s % wTop % array
zMid => s % zMid % array
h_edge => s % h_edge % array
+ viscosity => s % viscosity % array
vorticity => s % vorticity % array
divergence => s % divergence % array
ke => s % ke % array
@@ -237,7 +238,7 @@
! strictly only valid for config_h_mom_eddy_visc2 == constant
!
call mpas_timer_start("hmix", .false., velHmixTimer)
- call ocn_vel_hmix_tend(grid, divergence, vorticity, tend_u, err)
+ call ocn_vel_hmix_tend(grid, divergence, vorticity, viscosity, tend_u, err)
call mpas_timer_stop("hmix", velHmixTimer)
!
Modified: trunk/mpas/src/core_ocean/mpas_ocn_vel_hmix.F
===================================================================
--- trunk/mpas/src/core_ocean/mpas_ocn_vel_hmix.F 2012-10-22 19:14:22 UTC (rev 2241)
+++ trunk/mpas/src/core_ocean/mpas_ocn_vel_hmix.F 2012-10-22 19:36:32 UTC (rev 2242)
@@ -20,6 +20,7 @@
use mpas_configure
use mpas_timer
use ocn_vel_hmix_del2
+ use ocn_vel_hmix_leith
use ocn_vel_hmix_del4
implicit none
@@ -47,7 +48,7 @@
!
!--------------------------------------------------------------------
- type (timer_node), pointer :: del2Timer, del4Timer
+ type (timer_node), pointer :: del2Timer, leithTimer, del4Timer
!***********************************************************************
@@ -72,7 +73,7 @@
!
!-----------------------------------------------------------------------
- subroutine ocn_vel_hmix_tend(grid, divergence, vorticity, tend, err)!{{{
+ subroutine ocn_vel_hmix_tend(grid, divergence, vorticity, viscosity, tend, err)!{{{
!-----------------------------------------------------------------
!
@@ -98,6 +99,9 @@
real (kind=RKIND), dimension(:,:), intent(inout) :: &
tend !< Input/Output: velocity tendency
+ real (kind=RKIND), dimension(:,:), intent(inout) :: &
+ viscosity !< Input: viscosity
+
!-----------------------------------------------------------------
!
! output variables
@@ -112,7 +116,7 @@
!
!-----------------------------------------------------------------
- integer :: err1, err2
+ integer :: err1, err2, err3
!-----------------------------------------------------------------
!
@@ -122,14 +126,21 @@
!
!-----------------------------------------------------------------
+ viscosity = 0.0
+
call mpas_timer_start("del2", .false., del2Timer)
- call ocn_vel_hmix_del2_tend(grid, divergence, vorticity, tend, err1)
+ call ocn_vel_hmix_del2_tend(grid, divergence, vorticity, viscosity, tend, err1)
call mpas_timer_stop("del2", del2Timer)
+
+ call mpas_timer_start("leith", .false., leithTimer)
+ call ocn_vel_hmix_leith_tend(grid, divergence, vorticity, viscosity, tend, err2)
+ call mpas_timer_stop("leith", leithTimer)
+
call mpas_timer_start("del4", .false., del4Timer)
- call ocn_vel_hmix_del4_tend(grid, divergence, vorticity, tend, err2)
+ call ocn_vel_hmix_del4_tend(grid, divergence, vorticity, tend, err3)
call mpas_timer_stop("del4", del4Timer)
- err = ior(err1, err2)
+ err = ior(ior(err1, err2),err3)
!--------------------------------------------------------------------
@@ -163,12 +174,13 @@
integer, intent(out) :: err !< Output: error flag
- integer :: err1, err2
+ integer :: err1, err2, err3
call ocn_vel_hmix_del2_init(err1)
- call ocn_vel_hmix_del4_init(err2)
+ call ocn_vel_hmix_leith_init(err2)
+ call ocn_vel_hmix_del4_init(err3)
- err = ior(err1, err2)
+ err = ior(ior(err1, err2),err3)
!--------------------------------------------------------------------
Modified: trunk/mpas/src/core_ocean/mpas_ocn_vel_hmix_del2.F
===================================================================
--- trunk/mpas/src/core_ocean/mpas_ocn_vel_hmix_del2.F 2012-10-22 19:14:22 UTC (rev 2241)
+++ trunk/mpas/src/core_ocean/mpas_ocn_vel_hmix_del2.F 2012-10-22 19:36:32 UTC (rev 2242)
@@ -70,7 +70,7 @@
!
!-----------------------------------------------------------------------
- subroutine ocn_vel_hmix_del2_tend(grid, divergence, vorticity, tend, err)!{{{
+ subroutine ocn_vel_hmix_del2_tend(grid, divergence, vorticity, viscosity, tend, err)!{{{
!-----------------------------------------------------------------
!
@@ -96,6 +96,8 @@
real (kind=RKIND), dimension(:,:), intent(inout) :: &
tend !< Input/Output: velocity tendency
+ real (kind=RKIND), dimension(:,:), intent(inout) :: &
+ viscosity !< Input: viscosity
!-----------------------------------------------------------------
!
@@ -111,12 +113,11 @@
!
!-----------------------------------------------------------------
- integer :: iEdge, nEdgesSolve, cell1, cell2, vertex1, vertex2
- integer :: k
+ integer :: iEdge, nEdgesSolve, cell1, cell2, vertex1, vertex2, k
integer, dimension(:), pointer :: maxLevelEdgeTop
integer, dimension(:,:), pointer :: cellsOnEdge, verticesOnEdge, edgeMask
- real (kind=RKIND) :: u_diffusion, invLength1, invLength2
+ real (kind=RKIND) :: u_diffusion, invLength1, invLength2, visc2
real (kind=RKIND), dimension(:), pointer :: meshScalingDel2, &
dcEdge, dvEdge
@@ -158,10 +159,12 @@
-viscVortCoef &
*( vorticity(k,vertex2) - vorticity(k,vertex1) ) * invLength2
- u_diffusion = meshScalingDel2(iEdge) * eddyVisc2 * u_diffusion
+ visc2 = meshScalingDel2(iEdge) * eddyVisc2
- tend(k,iEdge) = tend(k,iEdge) + edgeMask(k, iEdge) * u_diffusion
+ tend(k,iEdge) = tend(k,iEdge) + edgeMask(k, iEdge) * visc2 * u_diffusion
+ viscosity(k,iEdge) = viscosity(k,iEdge) + visc2
+
end do
end do
Copied: trunk/mpas/src/core_ocean/mpas_ocn_vel_hmix_leith.F (from rev 2241, branches/ocean_projects/leith_mrp/src/core_ocean/mpas_ocn_vel_hmix_leith.F)
===================================================================
--- trunk/mpas/src/core_ocean/mpas_ocn_vel_hmix_leith.F (rev 0)
+++ trunk/mpas/src/core_ocean/mpas_ocn_vel_hmix_leith.F 2012-10-22 19:36:32 UTC (rev 2242)
@@ -0,0 +1,235 @@
+!|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
+!
+! ocn_vel_hmix_leith
+!
+!> \brief Ocean horizontal mixing - Leith parameterization
+!> \author Mark Petersen
+!> \date 22 October 2012
+!> \version SVN:$Id:$
+!> \details
+!> This module contains routines for computing horizontal mixing
+!> tendencies using the Leith parameterization.
+!
+!-----------------------------------------------------------------------
+
+module ocn_vel_hmix_leith
+
+ use mpas_grid_types
+ use mpas_configure
+
+ implicit none
+ private
+ save
+
+ !--------------------------------------------------------------------
+ !
+ ! Public parameters
+ !
+ !--------------------------------------------------------------------
+
+ !--------------------------------------------------------------------
+ !
+ ! Public member functions
+ !
+ !--------------------------------------------------------------------
+
+ public :: ocn_vel_hmix_leith_tend, &
+ ocn_vel_hmix_leith_init
+
+ !-------------------------------------------------------------------
+ !
+ ! Private module variables
+ !
+ !--------------------------------------------------------------------
+
+ logical :: hmixLeithOn !< integer flag to determine whether leith chosen
+
+ real (kind=RKIND) :: &
+ viscVortCoef
+
+!***********************************************************************
+
+contains
+
+!***********************************************************************
+!
+! routine ocn_vel_hmix_leith_tend
+!
+!> \brief Computes tendency term for horizontal momentum mixing with Leith parameterization
+!> \author Mark Petersen, Todd Ringler
+!> \date 22 October 2012
+!> \version SVN:$Id$
+!> \details
+!> This routine computes the horizontal mixing tendency for momentum
+!> based on the Leith closure. The Leith closure is the
+!> enstrophy-cascade analogy to the Smagorinsky (1963) energy-cascade
+!> closure, i.e. Leith (1996) assumes an inertial range of enstrophy flux
+!> moving toward the grid scale. The assumption of an enstrophy cascade
+!> and dimensional analysis produces right-hand-side dissipation,
+!> $\bf{D}$, of velocity of the form
+!> $ {\bf D} = </font>
<font color="black">abla \cdot \left( </font>
<font color="black">u_\ast </font>
<font color="blue">abla {\bf u} \right)
+!> = </font>
<font color="black">abla \cdot \left( \gamma \left| </font>
<font color="blue">abla \omega \right|
+!> \left( \Delta x \right)^3 </font>
<font color="blue">abla \bf{u} \right)
+!> where $\omega$ is the relative vorticity and $\gamma$ is a non-dimensional,
+!> $O(1)$ parameter. We set $\gamma=1$.
+
+!
+!-----------------------------------------------------------------------
+
+ subroutine ocn_vel_hmix_leith_tend(grid, divergence, vorticity, viscosity, 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
+
+ real (kind=RKIND), dimension(:,:), intent(inout) :: &
+ viscosity !< Input: viscosity
+
+ !-----------------------------------------------------------------
+ !
+ ! output variables
+ !
+ !-----------------------------------------------------------------
+
+ integer, intent(out) :: err !< Output: error flag
+
+ !-----------------------------------------------------------------
+ !
+ ! local variables
+ !
+ !-----------------------------------------------------------------
+
+ integer :: iEdge, nEdgesSolve, cell1, cell2, vertex1, vertex2, k
+ integer, dimension(:), pointer :: maxLevelEdgeTop
+ integer, dimension(:,:), pointer :: cellsOnEdge, verticesOnEdge, edgeMask
+
+ real (kind=RKIND) :: u_diffusion, invLength1, invLength2, visc2
+ real (kind=RKIND), dimension(:), pointer :: meshScaling, &
+ dcEdge, dvEdge
+
+ !-----------------------------------------------------------------
+ !
+ ! exit if this mixing is not selected
+ !
+ !-----------------------------------------------------------------
+
+ err = 0
+
+ if(.not.hmixLeithOn) return
+
+ nEdgesSolve = grid % nEdgesSolve
+ maxLevelEdgeTop => grid % maxLevelEdgeTop % array
+ cellsOnEdge => grid % cellsOnEdge % array
+ verticesOnEdge => grid % verticesOnEdge % array
+ meshScaling => grid % meshScaling % array
+ edgeMask => grid % edgeMask % 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)
+
+ invLength1 = 1.0 / dcEdge(iEdge)
+ invLength2 = 1.0 / dvEdge(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="blue">abla vorticity pointing from cell1 to cell2.
+
+ u_diffusion = ( divergence(k,cell2) - divergence(k,cell1) ) * invLength1 &
+ -viscVortCoef &
+ *( vorticity(k,vertex2) - vorticity(k,vertex1) ) * invLength2
+
+ ! Here the first line is (\delta x)^3
+ ! the second line is |</font>
<font color="blue">abla \omega|
+ ! and u_diffusion is </font>
<font color="blue">abla^2 u (see formula for $\bf{D}$ above).
+ visc2 = ( config_leith_parameter * config_leith_dx * meshScaling(iEdge) / 3.14)**3 &
+ * abs( vorticity(k,vertex2) - vorticity(k,vertex1) ) * invLength1 * sqrt(3.0)
+ visc2 = min(visc2, config_leith_visc2_max)
+
+ tend(k,iEdge) = tend(k,iEdge) + edgeMask(k, iEdge) * visc2 * u_diffusion
+
+ viscosity(k,iEdge) = viscosity(k,iEdge) + visc2
+
+ end do
+ end do
+
+ !--------------------------------------------------------------------
+
+ end subroutine ocn_vel_hmix_leith_tend!}}}
+
+!***********************************************************************
+!
+! routine ocn_vel_hmix_leith_init
+!
+!> \brief Initializes ocean momentum horizontal mixing with Leith parameterization
+!> \author Mark Petersen
+!> \date 22 October 2012
+!> \version SVN:$Id$
+!> \details
+!> This routine initializes a variety of quantities related to
+!> Leith parameterization for horizontal momentum mixing in the ocean.
+!
+!-----------------------------------------------------------------------
+
+ subroutine ocn_vel_hmix_leith_init(err)!{{{
+
+
+ integer, intent(out) :: err !< Output: error flag
+
+ !--------------------------------------------------------------------
+ !
+ ! set some local module variables based on input config choices
+ !
+ !--------------------------------------------------------------------
+
+ err = 0
+
+ hmixLeithOn = .false.
+
+ if (config_use_leith_del2) then
+ hmixLeithOn = .true.
+
+ if (config_visc_vorticity_term) then
+ viscVortCoef = config_visc_vorticity_visc2_scale
+ else
+ viscVortCoef = 0.0
+ endif
+
+ endif
+
+ !--------------------------------------------------------------------
+
+ end subroutine ocn_vel_hmix_leith_init!}}}
+
+!***********************************************************************
+
+end module ocn_vel_hmix_leith
+
+!|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
+! vim: foldmethod=marker
</font>
</pre>