<p><b>mpetersen@lanl.gov</b> 2012-10-25 14:41:41 -0600 (Thu, 25 Oct 2012)</p><p>branch commit, restart_reproducibility<br>
<br>
The following three changes were required to have bit-for-bit restartibility:<br>
<br>
1. Call diagnostics before implicit vertical mixing (IVM)<br>
<br>
Variables like rho, h_edge, and ke_edge must be updated at the end of<br>
a timestep but before IVM. Right now we just added an additional call<br>
to ocn_diagnostics_solve just before IVM. After kpp is added, we<br>
could improve performance by only computing the diagnostics required<br>
for that particular IVM scheme beforehand. Diagnostics must be<br>
recomputed after IVM as well, since IVM alters u and tracers.<br>
<br>
This item prevented bit restartibility for RK4 and split explicit.<br>
<br>
2. Compute wTop with the correct advection velocity<br>
<br>
There are two velocities for advection:<br>
u, used by momentum equation<br>
uTransport = u + uCorr + uBolus, used by h and tracer equation<br>
<br>
wTop in the vertical advection must be computed with the correct<br>
horizontal velocity for each conservation equation. In this revision,<br>
I changed wTop to accept u as an input variable. wTop is called with<br>
the correct advection velocity just before computing the tendancies.<br>
<br>
This item prevented bit restartibility for split explicit because wTop<br>
was always computed with uTransport, except upon restart when it used<br>
u.<br>
<br>
3. Add uBtr to restart file.<br>
<br>
In split explicit timestepping, the velocity u is not resplit at the<br>
beginning of each timestep. Instead, we keep uBtr from the previous<br>
timestep. Note that uBcl may now include a barotropic component,<br>
because the weights h have changed. That is OK, because the<br>
GBtrForcing variable subtracts out the barotropic component from the<br>
baroclinic.<br>
<br>
This item prevented bit restartibility for split explicit because uBtr<br>
was computed upon restart. In could be solved by freshly splitting u<br>
at the beginning of each timestep. However, Stage 1 is set up to take<br>
care of that with the G term, and the additional splitting is not<br>
needed. Instead, we chose to save uBtr to the restart file for bit<br>
restartability, since this is only the addition of a 2D field.<br>
</p><hr noshade><pre><font color="gray">Index: branches/ocean_projects/restart_reproducibility
===================================================================
--- branches/ocean_projects/restart_reproducibility 2012-10-25 20:24:32 UTC (rev 2268)
+++ branches/ocean_projects/restart_reproducibility 2012-10-25 20:41:41 UTC (rev 2269)
Property changes on: branches/ocean_projects/restart_reproducibility
___________________________________________________________________
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
## -22,3 +23,4 ##
/branches/omp_blocks/multiple_blocks:1803-2084
/branches/source_renaming:1082-1113
/branches/time_manager:924-962
+/trunk/mpas:2239-2242
\ No newline at end of property
Index: branches/ocean_projects/restart_reproducibility/src/core_ocean
===================================================================
--- branches/ocean_projects/restart_reproducibility/src/core_ocean 2012-10-25 20:24:32 UTC (rev 2268)
+++ branches/ocean_projects/restart_reproducibility/src/core_ocean 2012-10-25 20:41:41 UTC (rev 2269)
Property changes on: branches/ocean_projects/restart_reproducibility/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
## -24,3 +25,4 ##
/branches/omp_blocks/openmp_test/src/core_ocean_elements:2161-2201
/branches/source_renaming/src/core_ocean:1082-1113
/branches/time_manager/src/core_ocean:924-962
+/trunk/mpas/src/core_ocean:2239-2242
\ No newline at end of property
Modified: branches/ocean_projects/restart_reproducibility/src/core_ocean/mpas_ocn_mpas_core.F
===================================================================
--- branches/ocean_projects/restart_reproducibility/src/core_ocean/mpas_ocn_mpas_core.F 2012-10-25 20:24:32 UTC (rev 2268)
+++ branches/ocean_projects/restart_reproducibility/src/core_ocean/mpas_ocn_mpas_core.F 2012-10-25 20:41:41 UTC (rev 2269)
@@ -114,7 +114,7 @@
call ocn_init_h_zstar(domain)
endif
- call ocn_init_split_timestep(domain)
+ if (.not.config_do_restart) call ocn_init_split_timestep(domain)
write (0,'(a,a10)'), ' Vertical grid type is: ',config_vert_grid_type
@@ -265,8 +265,6 @@
= block % state % time_levs(1) % state % u % array(:,:) &
+ block % state % time_levs(1) % state % uBolusGM % array(:,:)
- call ocn_wtop(block % state % time_levs(1) % state,block % state % time_levs(1) % state, mesh)
-
call ocn_compute_mesh_scaling(mesh)
call mpas_rbf_interp_initialize(mesh)
Modified: branches/ocean_projects/restart_reproducibility/src/core_ocean/mpas_ocn_tendency.F
===================================================================
--- branches/ocean_projects/restart_reproducibility/src/core_ocean/mpas_ocn_tendency.F 2012-10-25 20:24:32 UTC (rev 2268)
+++ branches/ocean_projects/restart_reproducibility/src/core_ocean/mpas_ocn_tendency.F 2012-10-25 20:41:41 UTC (rev 2269)
@@ -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)
!
@@ -853,13 +854,43 @@
!> This routine computes the vertical velocity in the top layer for the ocean
!
!-----------------------------------------------------------------------
- subroutine ocn_wtop(s1,s2, grid)!{{{
- implicit none
+ subroutine ocn_wtop(grid,h,h_edge,u,wTop, err)!{{{
- type (state_type), intent(inout) :: s1 !< Input/Output: State 1 information
- type (state_type), intent(inout) :: s2 !< Input/Output: State 2 information
- type (mesh_type), intent(in) :: grid !< Input: Grid information
+ !-----------------------------------------------------------------
+ !
+ ! input variables
+ !
+ !-----------------------------------------------------------------
+ type (mesh_type), intent(in) :: &
+ grid !< Input: grid information
+
+ real (kind=RKIND), dimension(:,:), intent(in) :: &
+ h !< Input: thickness
+
+ real (kind=RKIND), dimension(:,:), intent(in) :: &
+ h_edge !< Input: h interpolated to an edge
+
+ real (kind=RKIND), dimension(:,:), intent(in) :: &
+ u !< Input: velocity
+
+ !-----------------------------------------------------------------
+ !
+ ! output variables
+ !
+ !-----------------------------------------------------------------
+
+ real (kind=RKIND), dimension(:,:), intent(out) :: &
+ wTop !< Output: vertical transport at top edge
+
+ integer, intent(out) :: err !< Output: error flag
+
+ !-----------------------------------------------------------------
+ !
+ ! local variables
+ !
+ !-----------------------------------------------------------------
+
integer :: iEdge, iCell, iVertex, k, cell1, cell2, vertex1, vertex2, eoe, i, j, cov
real (kind=RKIND) :: flux, vorticity_abs, h_vertex, workpv, rho0Inv, hSum, invAreaCell
@@ -868,7 +899,6 @@
real (kind=RKIND), dimension(:), pointer :: &
dvEdge, areaCell, zstarWeight
- real (kind=RKIND), dimension(:,:), pointer :: uTransport,h,wTop, h_edge
real (kind=RKIND), dimension(:), allocatable:: div_hu, h_tend_col
real (kind=RKIND) :: div_hu_btr
@@ -879,10 +909,7 @@
maxLevelCell, maxLevelEdgeTop, maxLevelEdgeBot, &
maxLevelVertexBot, maxLevelVertexTop
- h => s1 % h % array
- h_edge => s1 % h_edge % array
- uTransport => s2 % uTransport % array
- wTop => s2 % wTop % array
+ err = 0
nEdgesOnCell => grid % nEdgesOnCell % array
areaCell => grid % areaCell % array
@@ -922,7 +949,7 @@
iEdge = edgesOnCell(i, iCell)
do k = 1, maxLevelEdgeBot(iEdge)
- flux = uTransport(k, iEdge) * dvEdge(iEdge) * h_edge(k, iEdge)
+ flux = u(k, iEdge) * dvEdge(iEdge) * h_edge(k, iEdge)
flux = edgeSignOnCell(i, iCell) * flux * invAreaCell
div_hu(k) = div_hu(k) - flux
div_hu_btr = div_hu_btr - flux
Modified: branches/ocean_projects/restart_reproducibility/src/core_ocean/mpas_ocn_time_integration_rk4.F
===================================================================
--- branches/ocean_projects/restart_reproducibility/src/core_ocean/mpas_ocn_time_integration_rk4.F 2012-10-25 20:24:32 UTC (rev 2268)
+++ branches/ocean_projects/restart_reproducibility/src/core_ocean/mpas_ocn_time_integration_rk4.F 2012-10-25 20:41:41 UTC (rev 2269)
@@ -149,17 +149,19 @@
block => domain % blocklist
do while (associated(block))
- ! mrp 111206 put ocn_wtop call at top for ALE
- call ocn_wtop(block % provis, block % provis, block % mesh)
-
if (.not.config_implicit_vertical_mix) then
call ocn_vmix_coefs(block % mesh, block % provis, block % diagnostics, err)
end if
- call ocn_tend_h(block % tend, block % provis, block % mesh)
+
+ ! advection of u uses u, while advection of h and tracers use uTransport.
+ call ocn_wtop(block % mesh, block % provis % h % array, block % provis % h_edge % array, &
+ block % provis % u % array, block % provis % wTop % array, err)
call ocn_tend_u(block % tend, block % provis, block % diagnostics, block % mesh)
- ! mrp 110718 filter btr mode out of u_tend
- ! still got h perturbations with just this alone. Try to set uBtr=0 after full u computation
+ call ocn_wtop(block % mesh, block % provis % h % array, block % provis % h_edge % array, &
+ block % provis % uTransport % array, block % provis % wTop % array, err)
+ call ocn_tend_h(block % tend, block % provis, block % mesh)
+
if (config_rk_filter_btr_mode) then
call ocn_filter_btr_mode_tend_u(block % tend, block % provis, block % mesh)
endif
@@ -273,6 +275,15 @@
call mpas_timer_start("RK4-implicit vert mix")
block => domain % blocklist
do while(associated(block))
+
+ ! Call ocean diagnostic solve in preparation for vertical mixing. Note
+ ! it is called again after vertical mixing, because u and tracers change.
+ ! For Richardson vertical mixing, only rho, h_edge, and ke_edge need to
+ ! be computed. For kpp, more variables may be needed. Either way, this
+ ! could be made more efficient by only computing what is needed for the
+ ! implicit vmix routine that follows. mrp 121023.
+ call ocn_diagnostic_solve(dt, block % state % time_levs(2) % state, block % mesh)
+
call ocn_vmix_implicit(dt, block % mesh, block % diagnostics, block % state % time_levs(2) % state, err)
block => block % next
end do
Modified: branches/ocean_projects/restart_reproducibility/src/core_ocean/mpas_ocn_time_integration_split.F
===================================================================
--- branches/ocean_projects/restart_reproducibility/src/core_ocean/mpas_ocn_time_integration_split.F 2012-10-25 20:24:32 UTC (rev 2268)
+++ branches/ocean_projects/restart_reproducibility/src/core_ocean/mpas_ocn_time_integration_split.F 2012-10-25 20:41:41 UTC (rev 2269)
@@ -86,9 +86,9 @@
type (dm_info) :: dminfo
integer :: iCell, i,k,j, iEdge, cell1, cell2, split_explicit_step, split, &
eoe, oldBtrSubcycleTime, newBtrSubcycleTime, uPerpTime, BtrCorIter, &
- n_bcl_iter(config_n_ts_iter)
+ n_bcl_iter(config_n_ts_iter), stage1_tend_time
type (block_type), pointer :: block
- real (kind=RKIND) :: uhSum, hSum, flux, sshEdge, &
+ real (kind=RKIND) :: uhSum, hSum, flux, sshEdge, hEdge1, &
CoriolisTerm, uCorr, temp, temp_h, coef, FBtr_coeff, sshCell1, sshCell2
integer :: num_tracers, ucorr_coef, err
real (kind=RKIND), dimension(:,:), pointer :: &
@@ -117,6 +117,9 @@
! The baroclinic velocity needs be recomputed at the beginning of a
! timestep because the implicit vertical mixing is conducted on the
! total u. We keep uBtr from the previous timestep.
+ ! Note that uBcl may now include a barotropic component, because the
+ ! weights h have changed. That is OK, because the GBtrForcing variable
+ ! subtracts out the barotropic component from the baroclinic.
block % state % time_levs(1) % state % uBcl % array(k,iEdge) &
= block % state % time_levs(1) % state % u % array(k,iEdge) &
- block % state % time_levs(1) % state % uBtr % array( iEdge)
@@ -181,10 +184,22 @@
block => domain % blocklist
do while (associated(block))
+
+ stage1_tend_time = min(split_explicit_step,2)
+
if (.not.config_implicit_vertical_mix) then
- call ocn_vmix_coefs(block % mesh, block % state % time_levs(2) % state, block % diagnostics, err)
+ call ocn_vmix_coefs(block % mesh, block % state % time_levs(stage1_tend_time) % state, block % diagnostics, err)
end if
- call ocn_tend_u(block % tend, block % state % time_levs(2) % state, block % diagnostics, block % mesh)
+
+ ! compute wTop. Use u (rather than uTransport) for momentum advection.
+ ! Use the most recent time level available.
+ call ocn_wtop(block % mesh, block % state % time_levs(stage1_tend_time) % state % h % array, &
+ block % state % time_levs(stage1_tend_time) % state % h_edge % array, &
+ block % state % time_levs(stage1_tend_time) % state % u % array, &
+ block % state % time_levs(stage1_tend_time) % state % wTop % array, err)
+
+ call ocn_tend_u(block % tend, block % state % time_levs(stage1_tend_time) % state, block % diagnostics, block % mesh)
+
block => block % next
end do
@@ -246,6 +261,7 @@
= 0.5*( &
block % state % time_levs(1) % state % uBcl % array(k,iEdge) &
+ uTemp(k) - dt * block % state % time_levs(1) % state % GBtrForcing % array(iEdge))
+
enddo
enddo ! iEdge
@@ -767,8 +783,14 @@
! dwj: 02/22/12 splitting thickness and tracer tendency computations and halo updates to allow monotonic advection.
block => domain % blocklist
do while (associated(block))
- call ocn_wtop(block % state % time_levs(1) % state,block % state % time_levs(2) % state, block % mesh)
+ ! compute wTop. Use uTransport for advection of h and tracers.
+ ! Use time level 1 values of h and h_edge because h has not yet been computed for time level 2.
+ call ocn_wtop(block % mesh, block % state % time_levs(1) % state % h % array, &
+ block % state % time_levs(1) % state % h_edge % array, &
+ block % state % time_levs(2) % state % uTransport % array, &
+ block % state % time_levs(2) % state % wTop % array, err)
+
call ocn_tend_h(block % tend, block % state % time_levs(2) % state, block % mesh)
block => block % next
end do
@@ -921,7 +943,17 @@
call mpas_timer_start("se implicit vert mix")
block => domain % blocklist
do while(associated(block))
+
+ ! Call ocean diagnostic solve in preparation for vertical mixing. Note
+ ! it is called again after vertical mixing, because u and tracers change.
+ ! For Richardson vertical mixing, only rho, h_edge, and ke_edge need to
+ ! be computed. For kpp, more variables may be needed. Either way, this
+ ! could be made more efficient by only computing what is needed for the
+ ! implicit vmix routine that follows. mrp 121023.
+ call ocn_diagnostic_solve(dt, block % state % time_levs(2) % state, block % mesh)
+
call ocn_vmix_implicit(dt, block % mesh, block % diagnostics, block % state % time_levs(2) % state, err)
+
block => block % next
end do
@@ -959,6 +991,12 @@
= block % state % time_levs(2) % state % u % array(:,:) &
+ block % state % time_levs(2) % state % uBolusGM % array(:,:)
+!!!!!!!!!!!!!!!!!!!!!!!!!! mrp
+ ! Recompute wTop with all updated information. The wTop in stage 3 used
+ ! the previous h and u before implicit vertical mixing.
+!!!!!!
+! remove: call ocn_wtop(block % state % time_levs(2) % state,block % state % time_levs(2) % state, block % mesh)
+
call mpas_reconstruct(block % mesh, block % state % time_levs(2) % state % u % array, &
block % state % time_levs(2) % state % uReconstructX % array, &
block % state % time_levs(2) % state % uReconstructY % array, &
Copied: branches/ocean_projects/restart_reproducibility/src/core_ocean/mpas_ocn_vel_hmix_leith.F (from rev 2242, trunk/mpas/src/core_ocean/mpas_ocn_vel_hmix_leith.F)
===================================================================
--- branches/ocean_projects/restart_reproducibility/src/core_ocean/mpas_ocn_vel_hmix_leith.F (rev 0)
+++ branches/ocean_projects/restart_reproducibility/src/core_ocean/mpas_ocn_vel_hmix_leith.F 2012-10-25 20:41:41 UTC (rev 2269)
@@ -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>