<p><b>dwj07@fsu.edu</b> 2011-09-20 14:17:42 -0600 (Tue, 20 Sep 2011)</p><p><br>
-- BRANCH COMMIT --<br>
<br>
Adding modules for vertical mixing of momentum and tracers.<br>
<br>
These include options of explicit, and implicit mixing. As well<br>
as coefficients for constant, tanh, and richarson number mixing.<br>
<br>
In addition, calls to the new OcnEquationOfState module were commented out<br>
and reaplced with the older equation_of_state, as the new one is causing issues.<br>
<br>
This will hopefully be fixed in the next commit.<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-20 17:15:08 UTC (rev 1012)
+++ branches/ocean_projects/performance/src/core_ocean/Makefile 2011-09-20 20:17:42 UTC (rev 1013)
@@ -33,6 +33,10 @@
module_OcnEquationOfState.o \
module_OcnEquationOfStateLinear.o \
module_OcnEquationOfStateJM.o \
+ module_OcnVmix.o \
+ module_OcnVmixCoefsConst.o \
+ module_OcnVmixCoefsRich.o \
+ module_OcnVmixCoefsTanh.o \
module_time_integration.o \
module_global_diagnostics.o
@@ -109,6 +113,14 @@
module_OcnEquationOfStateJM.o:
+module_OcnVmix.o: module_OcnVmixCoefsConst.o module_OcnVmixCoefsRich.o module_OcnVmixCoefsTanh.o
+
+module_OcnVmixCoefsConst.o:
+
+module_OcnVmixCoefsRich.o:
+
+module_OcnVmixCoefsTanh.o:
+
module_mpas_core.o: module_advection.o \
module_OcnThickHadv.o \
module_OcnThickVadv.o \
@@ -142,6 +154,10 @@
module_OcnEquationOfState.o \
module_OcnEquationOfStateLinear.o \
module_OcnEquationOfStateJM.o \
+ module_OcnVmix.o \
+ module_OcnVmixCoefsConst.o \
+ module_OcnVmixCoefsRich.o \
+ module_OcnVmixCoefsTanh.o \
module_time_integration.o
clean:
Added: branches/ocean_projects/performance/src/core_ocean/module_OcnVmix.F
===================================================================
--- branches/ocean_projects/performance/src/core_ocean/module_OcnVmix.F (rev 0)
+++ branches/ocean_projects/performance/src/core_ocean/module_OcnVmix.F 2011-09-20 20:17:42 UTC (rev 1013)
@@ -0,0 +1,724 @@
+!|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
+!
+! OcnVmix
+!
+!> \brief MPAS ocean vertical mixing driver
+!> \author Doug Jacobsen
+!> \date 19 September 2011
+!> \version SVN:$Id:$
+!> \details
+!> This module is the main driver for
+!> vertical mixing in the ocean.
+!>
+!
+!-----------------------------------------------------------------------
+
+module OcnVmix
+
+ use grid_types
+ use configure
+ use timer
+
+ use OcnVmixCoefsConst
+ use OcnVmixCoefsTanh
+ use OcnVmixCoefsRich
+
+ implicit none
+ private
+ save
+
+ !--------------------------------------------------------------------
+ !
+ ! Public parameters
+ !
+ !--------------------------------------------------------------------
+
+ !--------------------------------------------------------------------
+ !
+ ! Public member functions
+ !
+ !--------------------------------------------------------------------
+
+ private :: tridiagonal_solve, &
+ tridiagonal_solve_mult
+
+ public :: OcnVmixCoefs, &
+ OcnVelVmixTendExplicit, &
+ OcnTracerVmixTendExplicit, &
+ OcnVelVmixTendImplicit, &
+ OcnTracerVmixTendImplicit, &
+ OcnVmixInit
+
+ !--------------------------------------------------------------------
+ !
+ ! Private module variables
+ !
+ !--------------------------------------------------------------------
+
+ logical :: explicitOn, implicitOn
+
+!***********************************************************************
+
+contains
+
+!***********************************************************************
+!
+! routine OcnVmixCoefs
+!
+!> \brief Computes coefficients for vertical mixing
+!> \author Doug Jacobsen
+!> \date 19 September 2011
+!> \version SVN:$Id$
+!> \details
+!> This routine computes the vertical mixing coefficients for momentum
+!> and tracers based user choices of mixing parameterization.
+!
+!-----------------------------------------------------------------------
+
+ subroutine OcnVmixCoefs(grid, s, d, err)!{{{
+
+ !-----------------------------------------------------------------
+ !
+ ! input variables
+ !
+ !-----------------------------------------------------------------
+
+ type (mesh_type), intent(in) :: &
+ grid !< Input: grid information
+
+ !-----------------------------------------------------------------
+ !
+ ! input/output variables
+ !
+ !-----------------------------------------------------------------
+
+ type (state_type), intent(inout) :: &
+ s !< Input/Output: state information
+
+ type (diagnostics_type), intent(inout) :: &
+ d !< Input/Output: diagnostic information
+
+ !-----------------------------------------------------------------
+ !
+ ! output variables
+ !
+ !-----------------------------------------------------------------
+
+ integer, intent(out) :: err
+
+ !-----------------------------------------------------------------
+ !
+ ! local variables
+ !
+ !-----------------------------------------------------------------
+
+ integer :: err1, err2, err3
+
+ !-----------------------------------------------------------------
+ !
+ ! call relevant routines for computing coefficients
+ !
+ !-----------------------------------------------------------------
+
+ err = 0
+
+ call OcnVmixCoefsConstBuild(grid, s, d, err1)
+ call OcnVmixCoefsTanhBuild(grid, s, d, err2)
+ call OcnVmixCoefsRichBuild(grid, s, d, err3)
+
+ err = err1 .or. err2 .or. err3
+
+ !--------------------------------------------------------------------
+
+ end subroutine OcnVmixCoefs!}}}
+
+!***********************************************************************
+!
+! routine OcnVelVmixTendExplict
+!
+!> \brief Computes tendencies for explict momentum vertical mixing
+!> \author Doug Jacobsen
+!> \date 19 September 2011
+!> \version SVN:$Id$
+!> \details
+!> This routine computes the tendencies for explicit vertical mixing for momentum
+!> using computed coefficients.
+!
+!-----------------------------------------------------------------------
+
+ subroutine OcnVelVmixTendExplicit(grid, u, h_edge, vertViscTopOfEdge, tend, err)!{{{
+
+ !-----------------------------------------------------------------
+ !
+ ! input variables
+ !
+ !-----------------------------------------------------------------
+
+ type (mesh_type), intent(in) :: &
+ grid !< Input: grid information
+
+ real (kind=RKIND), dimension(:,:), intent(in) :: &
+ u !< Input: velocity
+
+ real (kind=RKIND), dimension(:,:), intent(in) :: &
+ h_edge !< Input: thickness at edge
+
+ real (kind=RKIND), dimension(:,:), intent(in) :: &
+ vertViscTopOfEdge !< Input: vertical mixing coefficients
+
+ !-----------------------------------------------------------------
+ !
+ ! input/output variables
+ !
+ !-----------------------------------------------------------------
+
+ real (kind=RKIND), dimension(:,:), intent(inout) :: &
+ tend !< Input/Output: tendency information
+
+ !-----------------------------------------------------------------
+ !
+ ! output variables
+ !
+ !-----------------------------------------------------------------
+
+ integer, intent(out) :: err
+
+ !-----------------------------------------------------------------
+ !
+ ! local variables
+ !
+ !-----------------------------------------------------------------
+
+ integer :: iEdge, nEdgesSolve, k, nVertLevels
+
+ integer, dimension(:), pointer :: maxLevelEdgeTop
+
+ real (kind=RKIND), dimension(:), allocatable :: fluxVertTop
+
+ err = 0
+
+ if(implicitOn) return
+
+ call timer_start("compute_tend_u-explicit vert mix")
+
+ nEdgessolve = grid % nEdgesSolve
+ nVertLevels = grid % nVertLevels
+ maxLevelEdgeTop => grid % maxLevelEdgeTop % array
+
+ allocate(fluxVertTop(nVertLevels+1))
+ fluxVertTop(1) = 0.0
+ do iEdge=1,nEdgesSolve
+ do k=2,maxLevelEdgeTop(iEdge)
+ fluxVertTop(k) = vertViscTopOfEdge(k,iEdge) &
+ * ( u(k-1,iEdge) - u(k,iEdge) ) &
+ * 2 / (h_edge(k-1,iEdge) + h_edge(k,iEdge))
+ enddo
+ fluxVertTop(maxLevelEdgeTop(iEdge)+1) = 0.0
+
+ do k=1,maxLevelEdgeTop(iEdge)
+ tend(k,iEdge) = tend(k,iEdge) &
+ + (fluxVertTop(k) - fluxVertTop(k+1)) &
+ / h_edge(k,iEdge)
+ enddo
+
+ end do
+ deallocate(fluxVertTop)
+
+ call timer_stop("compute_tend_u-explicit vert mix")
+
+ !--------------------------------------------------------------------
+
+ end subroutine OcnVelVmixTendExplicit!}}}
+
+!***********************************************************************
+!
+! routine OcnVelVmixTendImplicit
+!
+!> \brief Computes tendencies for implicit momentum vertical mixing
+!> \author Doug Jacobsen
+!> \date 19 September 2011
+!> \version SVN:$Id$
+!> \details
+!> This routine computes the tendencies for implicit vertical mixing for momentum
+!> using computed coefficients.
+!
+!-----------------------------------------------------------------------
+
+ subroutine OcnVelVmixTendImplicit(grid, dt, ke_edge, vertViscTopOfEdge, h, h_edge, u, err)!{{{
+
+ !-----------------------------------------------------------------
+ !
+ ! input variables
+ !
+ !-----------------------------------------------------------------
+
+ type (mesh_type), intent(in) :: &
+ grid !< Input: grid information
+
+ real (kind=RKIND), dimension(:,:), intent(in) :: &
+ ke_edge !< Input: kinetic energy at edge
+
+ real (kind=RKIND), dimension(:,:), intent(in) :: &
+ vertViscTopOfEdge !< Input: vertical mixing coefficients
+
+ real (kind=RKIND), intent(in) :: &
+ dt !< Input: time step
+
+ real (kind=RKIND), dimension(:,:), intent(in) :: &
+ h !< Input: thickness at cell center
+
+ !-----------------------------------------------------------------
+ !
+ ! input/output variables
+ !
+ !-----------------------------------------------------------------
+
+ real (kind=RKIND), dimension(:,:), intent(inout) :: &
+ u !< Input: velocity
+
+ real (kind=RKIND), dimension(:,:), intent(inout) :: &
+ h_edge !< Input: thickness at edge
+
+ !-----------------------------------------------------------------
+ !
+ ! output variables
+ !
+ !-----------------------------------------------------------------
+
+ integer, intent(out) :: err
+
+ !-----------------------------------------------------------------
+ !
+ ! local variables
+ !
+ !-----------------------------------------------------------------
+
+ integer :: iEdge, nEdges, k, cell1, cell2, nVertLevels
+
+ integer, dimension(:), pointer :: maxLevelEdgeTop
+
+ integer, dimension(:,:), pointer :: cellsOnEdge
+
+ real (kind=RKIND), dimension(:), allocatable :: A, C, uTemp
+
+ err = 0
+
+ if(explicitOn) return
+
+ nEdges = grid % nEdges
+ nVertLevels = grid % nVertLevels
+ maxLevelEdgeTop => grid % maxLevelEdgeTop % array
+ cellsOnEdge => grid % cellsOnEdge % array
+
+ allocate(A(nVertLevels),C(nVertLevels),uTemp(nVertLevels))
+
+ do iEdge=1,nEdges
+ if (maxLevelEdgeTop(iEdge).gt.0) then
+
+ ! Compute A(k), C(k) for momentum
+ ! mrp 110315 efficiency note: for z-level, could precompute
+ ! -2.0*dt/(h(k)_h(k+1))/h(k) in setup
+ ! h_edge is computed in compute_solve_diag, and is not available yet.
+ ! This could be removed if hZLevel used instead.
+ cell1 = cellsOnEdge(1,iEdge)
+ cell2 = cellsOnEdge(2,iEdge)
+ do k=1,maxLevelEdgeTop(iEdge)
+ h_edge(k,iEdge) = 0.5 * (h(k,cell1) + h(k,cell2))
+ end do
+
+ do k=1,maxLevelEdgeTop(iEdge)-1
+ A(k) = -2.0*dt*vertViscTopOfEdge(k+1,iEdge) &
+ / (h_edge(k,iEdge) + h_edge(k+1,iEdge)) &
+ / h_edge(k,iEdge)
+ enddo
+ A(maxLevelEdgeTop(iEdge)) = -dt*config_bottom_drag_coeff &
+ *sqrt(2.0*ke_edge(k,iEdge))/h_edge(k,iEdge)
+
+ C(1) = 1 - A(1)
+ do k=2,maxLevelEdgeTop(iEdge)
+ C(k) = 1 - A(k) - A(k-1)
+ enddo
+
+ call tridiagonal_solve(A,C,A,u(:,iEdge),uTemp,maxLevelEdgeTop(iEdge))
+
+ u(1:maxLevelEdgeTop(iEdge),iEdge) = uTemp(1:maxLevelEdgeTop(iEdge))
+ u(maxLevelEdgeTop(iEdge)+1:nVertLevels,iEdge) = 0.0
+
+ end if
+ end do
+
+ !--------------------------------------------------------------------
+
+ end subroutine OcnVelVmixTendImplicit!}}}
+
+!***********************************************************************
+!
+! routine OcnTracerVmixTendExplict
+!
+!> \brief Computes tendencies for explict tracer vertical mixing
+!> \author Doug Jacobsen
+!> \date 19 September 2011
+!> \version SVN:$Id$
+!> \details
+!> This routine computes the tendencies for explicit vertical mixing for
+!> tracers using computed coefficients.
+!
+!-----------------------------------------------------------------------
+
+ subroutine OcnTracerVmixTendExplicit(grid, h, vertDiffTopOfCell, tracers, tend, err)!{{{
+
+ !-----------------------------------------------------------------
+ !
+ ! input variables
+ !
+ !-----------------------------------------------------------------
+
+ type (mesh_type), intent(in) :: &
+ grid !< Input: grid information
+
+ real (kind=RKIND), dimension(:,:), intent(in) :: &
+ h !< Input: thickness at cell center
+
+ real (kind=RKIND), dimension(:,:), intent(in) :: &
+ vertDiffTopOfCell !< Input: vertical mixing coefficients
+
+ real (kind=RKIND), dimension(:,:,:), intent(in) :: &
+ tracers !< Input: tracers
+
+ !-----------------------------------------------------------------
+ !
+ ! input/output variables
+ !
+ !-----------------------------------------------------------------
+
+ real (kind=RKIND), dimension(:,:,:), intent(inout) :: &
+ tend !< Input/Output: tendency information
+
+ !-----------------------------------------------------------------
+ !
+ ! output variables
+ !
+ !-----------------------------------------------------------------
+
+ integer, intent(out) :: err
+
+ !-----------------------------------------------------------------
+ !
+ ! local variables
+ !
+ !-----------------------------------------------------------------
+
+ integer :: iCell, nCellsSolve, k, iTracer, num_tracers, nVertLevels
+
+ integer, dimension(:), pointer :: maxLevelCell
+
+ real (kind=RKIND), dimension(:,:), allocatable :: fluxVertTop
+
+ err = 0
+
+ if(implicitOn) return
+
+ call timer_start("compute_scalar_tend-explicit vert diff")
+
+ nCellsSolve = grid % nCellsSolve
+ nVertLevels = grid % nVertLevels
+ num_tracers = size(tracers, dim=1)
+
+ maxLevelCell => grid % maxLevelCell % array
+
+ allocate(fluxVertTop(num_tracers,nVertLevels+1))
+ fluxVertTop(:,1) = 0.0
+ do iCell=1,nCellsSolve
+
+ do k=2,maxLevelCell(iCell)
+ do iTracer=1,num_tracers
+ ! compute \kappa_v d\phi/dz
+ fluxVertTop(iTracer,k) = vertDiffTopOfCell(k,iCell) &
+ * ( tracers(iTracer,k-1,iCell) &
+ - tracers(iTracer,k ,iCell) ) &
+ * 2 / (h(k-1,iCell) + h(k,iCell))
+
+ enddo
+ enddo
+ fluxVertTop(:,maxLevelCell(iCell)+1) = 0.0
+
+ do k=1,maxLevelCell(iCell)
+ do iTracer=1,num_tracers
+ ! This is h d/dz( fluxVertTop) but h and dz cancel, so
+ ! reduces to delta( fluxVertTop)
+ tend(iTracer,k,iCell) = tend(iTracer,k,iCell) &
+ + fluxVertTop(iTracer,k) - fluxVertTop(iTracer,k+1)
+ enddo
+ enddo
+!print '(a,50e12.2)', 'fluxVertTop',fluxVertTop(3,1:maxLevelCell(iCell)+1)
+!print '(a,50e12.2)', 'tend_tr ',tend_tr(3,1,1:maxLevelCell(iCell))
+ enddo ! iCell loop
+ deallocate(fluxVertTop)
+
+ call timer_stop("compute_scalar_tend-explicit vert diff")
+
+ !--------------------------------------------------------------------
+
+ end subroutine OcnTracerVmixTendExplicit!}}}
+
+!***********************************************************************
+!
+! routine OcnTracerVmixTendImplicit
+!
+!> \brief Computes tendencies for implicit tracer vertical mixing
+!> \author Doug Jacobsen
+!> \date 19 September 2011
+!> \version SVN:$Id$
+!> \details
+!> This routine computes the tendencies for implicit vertical mixing for
+!> tracers using computed coefficients.
+!
+!-----------------------------------------------------------------------
+
+ subroutine OcnTracerVmixTendImplicit(grid, dt, vertDiffTopOfCell, h, tracers, err)!{{{
+
+ !-----------------------------------------------------------------
+ !
+ ! input variables
+ !
+ !-----------------------------------------------------------------
+
+ type (mesh_type), intent(in) :: &
+ grid !< Input: grid information
+
+ real (kind=RKIND), dimension(:,:), intent(in) :: &
+ vertDiffTopOfCell !< Input: vertical mixing coefficients
+
+ real (kind=RKIND), intent(in) :: &
+ dt !< Input: time step
+
+ real (kind=RKIND), dimension(:,:), intent(in) :: &
+ h !< Input: thickness at cell center
+
+ !-----------------------------------------------------------------
+ !
+ ! input/output variables
+ !
+ !-----------------------------------------------------------------
+
+ real (kind=RKIND), dimension(:,:,:), intent(inout) :: &
+ tracers !< Input: tracers
+
+ !-----------------------------------------------------------------
+ !
+ ! output variables
+ !
+ !-----------------------------------------------------------------
+
+ integer, intent(out) :: err
+
+ !-----------------------------------------------------------------
+ !
+ ! local variables
+ !
+ !-----------------------------------------------------------------
+
+ integer :: iCell, nCells, k, nVertLevels, num_tracers
+
+ integer, dimension(:), pointer :: maxLevelCell
+
+ real (kind=RKIND), dimension(:), allocatable :: A, C
+ real (kind=RKIND), dimension(:,:), allocatable :: tracersTemp
+
+ err = 0
+
+ if(explicitOn) return
+
+ nCells = grid % nCells
+ nVertLevels = grid % nVertLevels
+ num_tracers = size(tracers, dim=1)
+ maxLevelCell => grid % maxLevelCell % array
+
+ allocate(A(nVertLevels),C(nVertLevels), tracersTemp(num_tracers,nVertLevels))
+
+ do iCell=1,nCells
+ ! Compute A(k), C(k) for tracers
+ ! mrp 110315 efficiency note: for z-level, could precompute
+ ! -2.0*dt/(h(k)_h(k+1))/h(k) in setup
+ do k=1,maxLevelCell(iCell)-1
+ A(k) = -2.0*dt*vertDiffTopOfCell(k+1,iCell) &
+ / (h(k,iCell) + h(k+1,iCell)) / h(k,iCell)
+ enddo
+
+ A(maxLevelCell(iCell)) = 0.0
+
+ C(1) = 1 - A(1)
+ do k=2,maxLevelCell(iCell)
+ C(k) = 1 - A(k) - A(k-1)
+ enddo
+
+ call tridiagonal_solve_mult(A,C,A,tracers(:,:,iCell), &
+ tracersTemp, maxLevelCell(iCell), nVertLevels,num_tracers)
+
+ tracers(:,1:maxLevelCell(iCell),iCell) = tracersTemp(:,1:maxLevelCell(iCell))
+ tracers(:,maxLevelCell(iCell)+1:nVertLevels,iCell) = -1e34
+ end do
+ deallocate(A,C,tracersTemp)
+
+
+ !--------------------------------------------------------------------
+
+ end subroutine OcnTracerVmixTendImplicit!}}}
+
+!***********************************************************************
+!
+! routine OcnVmixInit
+!
+!> \brief Initializes ocean vertical mixing quantities
+!> \author Doug Jacobsen
+!> \date 19 September 2011
+!> \version SVN:$Id$
+!> \details
+!> This routine initializes a variety of quantities related to
+!> vertical mixing in the ocean. This primarily determines if
+!> explicit or implicit vertical mixing is to be used.
+!
+!-----------------------------------------------------------------------
+
+
+ subroutine OcnVmixInit(err)!{{{
+
+ !--------------------------------------------------------------------
+
+ !-----------------------------------------------------------------
+ !
+ ! call individual init routines for each parameterization
+ !
+ !-----------------------------------------------------------------
+
+ integer, intent(out) :: err
+
+ integer :: err1, err2, err3
+
+ err = 0
+
+ explicitOn = .true.
+ implicitOn = .false.
+
+ if(config_implicit_vertical_mix) then
+ explicitOn = .false.
+ implicitOn =.true.
+ end if
+
+ call OcnVmixCoefsConstInit(err1)
+ call OcnVmixCoefsTanhInit(err2)
+ call OcnVmixCoefsRichInit(err3)
+
+ err = err .or. err1 .or. err2 .or. err3
+
+ !--------------------------------------------------------------------
+
+ end subroutine OcnVmixInit!}}}
+
+subroutine tridiagonal_solve(a,b,c,r,x,n)!{{{
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+! Solve the matrix equation Ax=r for x, where A is tridiagonal.
+! A is an nxn matrix, with:
+! a sub-diagonal, filled from 1:n-1 (a(1) appears on row 2)
+! b diagonal, filled from 1:n
+! c sup-diagonal, filled from 1:n-1 (c(1) apears on row 1)
+!
+! Input: a,b,c,r,n
+!
+! Output: x
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+ implicit none
+
+ integer,intent(in) :: n
+ real (KIND=RKIND), dimension(n), intent(in) :: a,b,c,r
+ real (KIND=RKIND), dimension(n), intent(out) :: x
+ real (KIND=RKIND), dimension(n) :: bTemp,rTemp
+ real (KIND=RKIND) :: m
+ integer i
+
+ call timer_start("tridiagonal_solve")
+
+ ! Use work variables for b and r
+ bTemp(1) = b(1)
+ rTemp(1) = r(1)
+
+ ! First pass: set the coefficients
+ do i = 2,n
+ m = a(i-1)/bTemp(i-1)
+ bTemp(i) = b(i) - m*c(i-1)
+ rTemp(i) = r(i) - m*rTemp(i-1)
+ end do
+
+ x(n) = rTemp(n)/bTemp(n)
+ ! Second pass: back-substition
+ do i = n-1, 1, -1
+ x(i) = (rTemp(i) - c(i)*x(i+1))/bTemp(i)
+ end do
+
+ call timer_stop("tridiagonal_solve")
+
+end subroutine tridiagonal_solve!}}}
+
+subroutine tridiagonal_solve_mult(a,b,c,r,x,n,nDim,nSystems)!{{{
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+! Solve the matrix equation Ax=r for x, where A is tridiagonal.
+! A is an nxn matrix, with:
+! a sub-diagonal, filled from 1:n-1 (a(1) appears on row 2)
+! b diagonal, filled from 1:n
+! c sup-diagonal, filled from 1:n-1 (c(1) apears on row 1)
+!
+! Input: a,b,c,r,n
+!
+! Output: x
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+ implicit none
+
+ integer,intent(in) :: n, nDim, nSystems
+ real (KIND=RKIND), dimension(n), intent(in) :: a,b,c
+ real (KIND=RKIND), dimension(nSystems,nDim), intent(in) :: r
+ real (KIND=RKIND), dimension(nSystems,nDim), intent(out) :: x
+ real (KIND=RKIND), dimension(n) :: bTemp
+ real (KIND=RKIND), dimension(nSystems,n) :: rTemp
+ real (KIND=RKIND) :: m
+ integer i,j
+
+ call timer_start("tridiagonal_solve_mult")
+
+ ! Use work variables for b and r
+ bTemp(1) = b(1)
+ do j = 1,nSystems
+ rTemp(j,1) = r(j,1)
+ end do
+
+ ! First pass: set the coefficients
+ do i = 2,n
+ m = a(i-1)/bTemp(i-1)
+ bTemp(i) = b(i) - m*c(i-1)
+ do j = 1,nSystems
+ rTemp(j,i) = r(j,i) - m*rTemp(j,i-1)
+ end do
+ end do
+
+ do j = 1,nSystems
+ x(j,n) = rTemp(j,n)/bTemp(n)
+ end do
+ ! Second pass: back-substition
+ do i = n-1, 1, -1
+ do j = 1,nSystems
+ x(j,i) = (rTemp(j,i) - c(i)*x(j,i+1))/bTemp(i)
+ end do
+ end do
+
+ call timer_stop("tridiagonal_solve_mult")
+
+end subroutine tridiagonal_solve_mult!}}}
+
+!***********************************************************************
+
+end module OcnVmix
+
+!|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
+
+!vim: foldmethod=marker
Added: branches/ocean_projects/performance/src/core_ocean/module_OcnVmixCoefsConst.F
===================================================================
--- branches/ocean_projects/performance/src/core_ocean/module_OcnVmixCoefsConst.F (rev 0)
+++ branches/ocean_projects/performance/src/core_ocean/module_OcnVmixCoefsConst.F 2011-09-20 20:17:42 UTC (rev 1013)
@@ -0,0 +1,306 @@
+!|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
+!
+! OcnVmixCoefsConst
+!
+!> \brief MPAS ocean vertical mixing coefficients
+!> \author Doug Jacobsen
+!> \date 19 September 2011
+!> \version SVN:$Id:$
+!> \details
+!> This module contains the routines for computing
+!> constant vertical mixing coefficients.
+!>
+!
+!-----------------------------------------------------------------------
+
+module OcnVmixCoefsConst
+
+ use grid_types
+ use configure
+ use timer
+
+ implicit none
+ private
+ save
+
+ !--------------------------------------------------------------------
+ !
+ ! Public parameters
+ !
+ !--------------------------------------------------------------------
+
+ !--------------------------------------------------------------------
+ !
+ ! Public member functions
+ !
+ !--------------------------------------------------------------------
+
+ private :: OcnVelVmixCoefsConst, &
+ OcnTracerVmixCoefsConst
+
+ public :: OcnVmixCoefsConstBuild, &
+ OcnVmixCoefsConstInit
+
+ !--------------------------------------------------------------------
+ !
+ ! Private module variables
+ !
+ !--------------------------------------------------------------------
+
+ logical :: constViscOn, constDiffOn
+
+ real (kind=RKIND) :: constVisc, constDiff
+
+
+!***********************************************************************
+
+contains
+
+!***********************************************************************
+!
+! routine OcnVmixCoefsConstBuild
+!
+!> \brief Computes coefficients for vertical mixing
+!> \author Doug Jacobsen
+!> \date 19 September 2011
+!> \version SVN:$Id$
+!> \details
+!> This routine computes the vertical mixing coefficients for momentum
+!> and tracers based user choices of mixing parameterization.
+!
+!-----------------------------------------------------------------------
+
+ subroutine OcnVmixCoefsConstBuild(grid, s, d, err)!{{{
+
+ !-----------------------------------------------------------------
+ !
+ ! input variables
+ !
+ !-----------------------------------------------------------------
+
+ type (mesh_type), intent(in) :: &
+ grid !< Input: grid information
+
+ !-----------------------------------------------------------------
+ !
+ ! input/output variables
+ !
+ !-----------------------------------------------------------------
+
+ type (state_type), intent(inout) :: &
+ s !< Input/Output: state information
+
+ type (diagnostics_type), intent(inout) :: &
+ d !< Input/Output: diagnostic information
+
+ !-----------------------------------------------------------------
+ !
+ ! output variables
+ !
+ !-----------------------------------------------------------------
+
+ integer, intent(out) :: err
+
+ !-----------------------------------------------------------------
+ !
+ ! local variables
+ !
+ !-----------------------------------------------------------------
+
+ integer :: err1, err2
+
+ real (kind=RKIND), dimension(:,:), pointer :: &
+ vertViscTopOfEdge, vertDiffTopOfCell
+
+ !-----------------------------------------------------------------
+ !
+ ! call relevant routines for computing tendencies
+ ! note that the user can choose multiple options and the
+ ! tendencies will be added together
+ !
+ !-----------------------------------------------------------------
+
+ err = 0
+ if((.not.constViscOn) .and. (.not.constDiffOn)) return
+
+ vertViscTopOfEdge => d % vertViscTopOfEdge % array
+ vertDiffTopOfCell => d % vertDiffTopOfCell % array
+
+ call OcnVelVmixCoefsConst(grid, vertViscTopOfEdge, err1)
+ call OcnTracerVmixCoefsConst(grid, vertDiffTopOfCell, err2)
+
+ err = err1 .or. err2
+
+ !--------------------------------------------------------------------
+
+ end subroutine OcnVmixCoefsConstBuild!}}}
+
+!***********************************************************************
+!
+! routine OcnVelVmixCoefsConst
+!
+!> \brief Computes coefficients for vertical momentum mixing
+!> \author Doug Jacobsen
+!> \date 19 September 2011
+!> \version SVN:$Id$
+!> \details
+!> This routine computes the constant vertical mixing coefficients for momentum
+!
+!-----------------------------------------------------------------------
+
+ subroutine OcnVelVmixCoefsConst(grid, vertViscTopOfEdge, err)!{{{
+
+ !-----------------------------------------------------------------
+ !
+ ! input variables
+ !
+ !-----------------------------------------------------------------
+
+ type (mesh_type), intent(in) :: &
+ grid !< Input: grid information
+
+ !-----------------------------------------------------------------
+ !
+ ! input/output variables
+ !
+ !-----------------------------------------------------------------
+
+ !-----------------------------------------------------------------
+ !
+ ! output variables
+ !
+ !-----------------------------------------------------------------
+
+ real (kind=RKIND), dimension(:,:), intent(out) :: vertViscTopOfEdge
+
+ integer, intent(out) :: err
+
+ !-----------------------------------------------------------------
+ !
+ ! local variables
+ !
+ !-----------------------------------------------------------------
+
+ err = 0
+
+ if(.not.constViscOn) return
+
+ vertViscTopOfEdge = constVisc
+
+ !--------------------------------------------------------------------
+
+ end subroutine OcnVelVmixCoefsConst!}}}
+
+!***********************************************************************
+!
+! routine OcnTracerVmixCoefsConst
+!
+!> \brief Computes coefficients for vertical tracer mixing
+!> \author Doug Jacobsen
+!> \date 19 September 2011
+!> \version SVN:$Id$
+!> \details
+!> This routine computes the constant vertical mixing coefficients for tracers
+!
+!-----------------------------------------------------------------------
+
+ subroutine OcnTracerVmixCoefsConst(grid, vertDiffTopOfCell, err)!{{{
+
+ !-----------------------------------------------------------------
+ !
+ ! input variables
+ !
+ !-----------------------------------------------------------------
+
+ type (mesh_type), intent(in) :: &
+ grid !< Input: grid information
+
+ !-----------------------------------------------------------------
+ !
+ ! input/output variables
+ !
+ !-----------------------------------------------------------------
+
+ !-----------------------------------------------------------------
+ !
+ ! output variables
+ !
+ !-----------------------------------------------------------------
+
+ real (kind=RKIND), dimension(:,:), intent(out) :: vertDiffTopOfCell
+
+ integer, intent(out) :: err
+
+ !-----------------------------------------------------------------
+ !
+ ! local variables
+ !
+ !-----------------------------------------------------------------
+
+ err = 0
+
+ if(.not.constDiffOn) return
+
+ vertDiffTopOfCell = constDiff
+
+ !--------------------------------------------------------------------
+
+ end subroutine OcnTracerVmixCoefsConst!}}}
+
+
+!***********************************************************************
+!
+! routine OcnVmixCoefsConstInit
+!
+!> \brief Initializes ocean momentum vertical mixing quantities
+!> \author Doug Jacobsen
+!> \date 19 September 2011
+!> \version SVN:$Id$
+!> \details
+!> This routine initializes a variety of quantities related to
+!> vertical velocity mixing in the ocean. Since a variety of
+!> parameterizations are available, this routine primarily calls the
+!> individual init routines for each parameterization.
+!
+!-----------------------------------------------------------------------
+
+
+ subroutine OcnVmixCoefsConstInit(err)!{{{
+
+ !--------------------------------------------------------------------
+
+ !-----------------------------------------------------------------
+ !
+ ! call individual init routines for each parameterization
+ !
+ !-----------------------------------------------------------------
+
+ integer, intent(out) :: err
+
+ err = 0
+
+ constViscOn = .false.
+ constDiffOn = .false.
+
+ if (config_vert_visc_type.eq.'const') then
+ constViscOn = .true.
+ constVisc = config_vert_visc
+ endif
+
+ if (config_vert_diff_type.eq.'const') then
+ constDiffOn = .true.
+ constDiff = config_vert_diff
+ endif
+
+
+ !--------------------------------------------------------------------
+
+ end subroutine OcnVmixCoefsConstInit!}}}
+
+!***********************************************************************
+
+end module OcnVmixCoefsConst
+
+!|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
+
+!vim: foldmethod=marker
Added: branches/ocean_projects/performance/src/core_ocean/module_OcnVmixCoefsRich.F
===================================================================
--- branches/ocean_projects/performance/src/core_ocean/module_OcnVmixCoefsRich.F (rev 0)
+++ branches/ocean_projects/performance/src/core_ocean/module_OcnVmixCoefsRich.F 2011-09-20 20:17:42 UTC (rev 1013)
@@ -0,0 +1,915 @@
+!|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
+!
+! OcnVmixCoefsRich
+!
+!> \brief MPAS ocean vertical mixing coefficients
+!> \author Doug Jacobsen
+!> \date 19 September 2011
+!> \version SVN:$Id:$
+!> \details
+!> This module contains the routines for computing
+!> richardson vertical mixing coefficients.
+!>
+!
+!-----------------------------------------------------------------------
+
+module OcnVmixCoefsRich
+
+ use grid_types
+ use configure
+ use constants
+ use timer
+
+ use OcnEquationOfState
+
+ implicit none
+ private
+ save
+
+ !--------------------------------------------------------------------
+ !
+ ! Public parameters
+ !
+ !--------------------------------------------------------------------
+
+ !--------------------------------------------------------------------
+ !
+ ! Public member functions
+ !
+ !--------------------------------------------------------------------
+
+ private :: OcnVelVmixCoefsRich, &
+ OcnTracerVmixCoefsRich, &
+ OcnVmixGetRichNumbers
+
+ public :: OcnVmixCoefsRichBuild, &
+ OcnVmixCoefsRichInit
+
+ !--------------------------------------------------------------------
+ !
+ ! Private module variables
+ !
+ !--------------------------------------------------------------------
+
+ logical :: richViscOn, richDiffOn
+
+!***********************************************************************
+
+contains
+
+!***********************************************************************
+!
+! routine OcnVmixCoefsRichBuild
+!
+!> \brief Computes coefficients for vertical mixing
+!> \author Doug Jacobsen
+!> \date 19 September 2011
+!> \version SVN:$Id$
+!> \details
+!> This routine computes the vertical mixing coefficients for momentum
+!> and tracers based user choices of mixing parameterization.
+!
+!-----------------------------------------------------------------------
+
+ subroutine OcnVmixCoefsRichBuild(grid, s, d, err)!{{{
+
+ !-----------------------------------------------------------------
+ !
+ ! input variables
+ !
+ !-----------------------------------------------------------------
+
+ type (mesh_type), intent(in) :: &
+ grid !< Input: grid information
+
+ !-----------------------------------------------------------------
+ !
+ ! input/output variables
+ !
+ !-----------------------------------------------------------------
+
+ type (state_type), intent(inout) :: &
+ s !< Input/Output: state information
+
+ type (diagnostics_type), intent(inout) :: &
+ d !< Input/Output: diagnostic information
+
+ !-----------------------------------------------------------------
+ !
+ ! output variables
+ !
+ !-----------------------------------------------------------------
+
+ integer, intent(out) :: err
+
+ !-----------------------------------------------------------------
+ !
+ ! local variables
+ !
+ !-----------------------------------------------------------------
+
+ integer :: err1, err2, err3, indexT, indexS
+
+ real (kind=RKIND), dimension(:,:), pointer :: &
+ vertViscTopOfEdge, vertDiffTopOfCell, u, h, h_edge, rho, rhoDisplaced
+
+ real (kind=RKIND), dimension(:,:), pointer :: RiTopOfEdge, RiTopOfCell
+
+ real (kind=RKIND), dimension(:,:,:), pointer :: tracers
+
+ !-----------------------------------------------------------------
+ !
+ ! call relevant routines for computing tendencies
+ ! note that the user can choose multiple options and the
+ ! tendencies will be added together
+ !
+ !-----------------------------------------------------------------
+
+ err = 0
+ if((.not.richViscOn) .and. (.not.richDiffOn)) return
+
+ indexT = s%index_temperature
+ indexS = s%index_salinity
+
+ vertViscTopOfEdge => d % vertViscTopOfEdge % array
+ vertDiffTopOfCell => d % vertDiffTopOfCell % array
+ RiTopOfEdge => d % RiTopOfEdge % array
+ RiTopOfCell => d % RiTopOfCell % array
+
+ u => s % u % array
+ h => s % h % array
+ h_edge => s % h_edge % array
+ rho => s % rho % array
+ rhoDisplaced => s % rhoDisplaced % array
+ tracers => s % tracers % array
+
+ call equation_of_state(s, grid, 0, 'relative')
+ call equation_of_state(s, grid, 1, 'relative')
+
+ call OcnVmixGetRichNumbers(grid, indexT, indexS, u, h, h_edge, &
+ rho, rhoDisplaced, tracers, RiTopOfEdge, RiTopOfCell, err1)
+
+ call OcnVelVmixCoefsRich(grid, RiTopOfEdge, h_edge, vertViscTopOfEdge, err2)
+ call OcnTracerVmixCoefsRich(grid, RiTopOfCell, h, vertDiffTopOfCell, err3)
+
+ err = err1 .or. err2 .or. err3
+
+ !--------------------------------------------------------------------
+
+ end subroutine OcnVmixCoefsRichBuild!}}}
+
+!***********************************************************************
+!
+! routine OcnVelVmixCoefsRich
+!
+!> \brief Computes coefficients for vertical momentum mixing
+!> \author Doug Jacobsen
+!> \date 19 September 2011
+!> \version SVN:$Id$
+!> \details
+!> This routine computes the richardson vertical mixing coefficients for momentum
+!
+!-----------------------------------------------------------------------
+
+ subroutine OcnVelVmixCoefsRich(grid, RiTopOfEdge, h_edge, vertViscTopOfEdge, err)!{{{
+
+ !-----------------------------------------------------------------
+ !
+ ! input variables
+ !
+ !-----------------------------------------------------------------
+
+ type (mesh_type), intent(in) :: &
+ grid !< Input: grid information
+
+ real (kind=RKIND), dimension(:,:), intent(in) :: &
+ h_edge !< Input: thickness at edge
+
+ real (kind=RKIND), dimension(:,:), intent(in) :: &
+ RiTopOfEdge !< Richardson number at top of edge
+
+ !-----------------------------------------------------------------
+ !
+ ! input/output variables
+ !
+ !-----------------------------------------------------------------
+
+ !-----------------------------------------------------------------
+ !
+ ! output variables
+ !
+ !-----------------------------------------------------------------
+
+ real (kind=RKIND), dimension(:,:), intent(out) :: vertViscTopOfEdge
+
+ integer, intent(out) :: err
+
+ !-----------------------------------------------------------------
+ !
+ ! local variables
+ !
+ !-----------------------------------------------------------------
+
+ integer :: iEdge, nEdges, k
+
+ integer, dimension(:), pointer :: maxLevelEdgeTop
+
+ err = 0
+
+ if(.not.richViscOn) return
+
+ nEdges = grid % nEdges
+
+ maxLevelEdgeTop => grid % maxLevelEdgeTop % array
+
+ vertViscTopOfEdge = 0.0
+ do iEdge = 1,nEdges
+ do k = 2,maxLevelEdgeTop(iEdge)
+ ! mrp 110324 efficiency note: this if is inside iEdge and k loops.
+ ! Perhaps there is a more efficient way to do this.
+ if (RiTopOfEdge(k,iEdge)>0.0) then
+ vertViscTopOfEdge(k,iEdge) = config_bkrd_vert_visc &
+ + config_rich_mix / (1.0 + 5.0*RiTopOfEdge(k,iEdge))**2
+ ! maltrud do limiting of coefficient--should not be necessary
+ ! also probably better logic could be found
+ if (vertViscTopOfEdge(k,iEdge) > config_convective_visc) then
+ if( config_implicit_vertical_mix) then
+ vertViscTopOfEdge(k,iEdge) = config_convective_visc
+ else
+ vertViscTopOfEdge(k,iEdge) = &
+ ((h_edge(k-1,iEdge)+h_edge(k,iEdge))/2.0)**2/config_dt/4.0
+ end if
+ end if
+ else
+ ! mrp 110324 efficiency note: this if is inside iCell and k loops.
+ if (config_implicit_vertical_mix) then
+ ! for Ri<0 and implicit mix, use convective diffusion
+ vertViscTopOfEdge(k,iEdge) = config_convective_visc
+ else
+ ! for Ri<0 and explicit vertical mix,
+ ! use maximum diffusion allowed by CFL criterion
+ ! mrp 110324 efficiency note: for z-level, could use fixed
+ ! grid array hMeanTopZLevel and compute maxdiff on startup.
+ vertViscTopOfEdge(k,iEdge) = &
+ ((h_edge(k-1,iEdge)+h_edge(k,iEdge))/2.0)**2/config_dt/4.0
+ end if
+ end if
+ end do
+ end do
+
+
+ !--------------------------------------------------------------------
+
+ end subroutine OcnVelVmixCoefsRich!}}}
+
+!***********************************************************************
+!
+! routine OcnTracerVmixCoefsRich
+!
+!> \brief Computes coefficients for vertical tracer mixing
+!> \author Doug Jacobsen
+!> \date 19 September 2011
+!> \version SVN:$Id$
+!> \details
+!> This routine computes the richardson vertical mixing coefficients for tracers
+!
+!-----------------------------------------------------------------------
+
+ subroutine OcnTracerVmixCoefsRich(grid, RiTopOfCell, h, vertDiffTopOfCell, err)!{{{
+
+ !-----------------------------------------------------------------
+ !
+ ! input variables
+ !
+ !-----------------------------------------------------------------
+
+ type (mesh_type), intent(in) :: &
+ grid !< Input: grid information
+
+ real (kind=RKIND), dimension(:,:), intent(in) :: &
+ h !< Input: thickness at cell center
+
+ real (kind=RKIND), dimension(:,:), intent(in) :: &
+ RiTopOfCell !< Input: Richardson number at top of cell
+
+ !-----------------------------------------------------------------
+ !
+ ! input/output variables
+ !
+ !-----------------------------------------------------------------
+
+ !-----------------------------------------------------------------
+ !
+ ! output variables
+ !
+ !-----------------------------------------------------------------
+
+ real (kind=RKIND), dimension(:,:), intent(out) :: vertDiffTopOfCell
+
+ integer, intent(out) :: err
+
+ !-----------------------------------------------------------------
+ !
+ ! local variables
+ !
+ !-----------------------------------------------------------------
+
+ integer :: iCell, nCells, k
+
+ integer, dimension(:), pointer :: maxLevelCell
+
+ real (kind=RKIND) :: coef
+
+ err = 0
+
+ if(.not.richDiffOn) return
+
+ nCells = grid % nCells
+
+ maxLevelCell => grid % maxLevelCell % array
+
+ vertDiffTopOfCell = 0.0
+ coef = -gravity/1000.0/2.0
+ do iCell = 1,nCells
+ do k = 2,maxLevelCell(iCell)
+ ! mrp 110324 efficiency note: this if is inside iCell and k loops.
+ ! Perhaps there is a more efficient way to do this.
+ if (RiTopOfCell(k,iCell)>0.0) then
+ vertDiffTopOfCell(k,iCell) = config_bkrd_vert_diff &
+ + (config_bkrd_vert_visc &
+ + config_rich_mix / (1.0 + 5.0*RiTopOfCell(k,iCell))**2) &
+ / (1.0 + 5.0*RiTopOfCell(k,iCell))
+ ! maltrud do limiting of coefficient--should not be necessary
+ ! also probably better logic could be found
+ if (vertDiffTopOfCell(k,iCell) > config_convective_diff) then
+ if (config_implicit_vertical_mix) then
+ vertDiffTopOfCell(k,iCell) = config_convective_diff
+ else
+ vertDiffTopOfCell(k,iCell) = &
+ ((h(k-1,iCell)+h(k,iCell))/2.0)**2/config_dt/4.0
+ end if
+ end if
+ else
+ ! mrp 110324 efficiency note: this if is inside iCell and k loops.
+ if (config_implicit_vertical_mix) then
+ ! for Ri<0 and implicit mix, use convective diffusion
+ vertDiffTopOfCell(k,iCell) = config_convective_diff
+ else
+ ! for Ri<0 and explicit vertical mix,
+ ! use maximum diffusion allowed by CFL criterion
+ ! mrp 110324 efficiency note: for z-level, could use fixed
+ ! grid array hMeanTopZLevel and compute maxdiff on startup.
+ vertDiffTopOfCell(k,iCell) = &
+ ((h(k-1,iCell)+h(k,iCell))/2.0)**2/config_dt/4.0
+ end if
+ end if
+ end do
+ end do
+
+
+ !--------------------------------------------------------------------
+
+ end subroutine OcnTracerVmixCoefsRich!}}}
+
+!***********************************************************************
+!
+! routine OcnVmixGetRichNumbers
+!
+!> \brief Build richardson numbers for vertical mixing
+!> \author Doug Jacobsen
+!> \date 19 September 2011
+!> \version SVN:$Id$
+!> \details
+!> This routine builds the arrays needed for richardson number vertical
+!> mixing coefficients.
+!
+!-----------------------------------------------------------------------
+
+ subroutine OcnVmixGetRichNumbers(grid, indexT, indexS, u, h, h_edge, & !{{{
+ rho, rhoDisplaced, tracers, RiTopOfEdge, RiTopOfCell, err)
+
+ !-----------------------------------------------------------------
+ !
+ ! input variables
+ !
+ !-----------------------------------------------------------------
+
+ type (mesh_type), intent(in) :: &
+ grid !< Input: grid information
+
+ integer, intent(in) :: indexT, indexS
+
+ real (kind=RKIND), dimension(:,:), intent(in) :: u, h, h_edge
+
+ real (kind=RKIND), dimension(:,:,:), intent(in) :: tracers
+
+ !-----------------------------------------------------------------
+ !
+ ! input/output variables
+ !
+ !-----------------------------------------------------------------
+
+ !-----------------------------------------------------------------
+ !
+ ! output variables
+ !
+ !-----------------------------------------------------------------
+
+ real (kind=RKIND), dimension(:,:), intent(inout) :: rho, rhoDisplaced, &
+ RiTopOfEdge, RiTopOfCell
+
+ integer, intent(inout) :: err
+
+ !-----------------------------------------------------------------
+ !
+ ! local variables
+ !
+ !-----------------------------------------------------------------
+
+ integer :: nVertLevels, nCells, nEdges, iCell, iEdge, k
+ integer :: cell1, cell2
+
+ integer, dimension(:), pointer :: maxLevelCell, maxLevelEdgeTop, maxLevelEdgeBot
+ integer, dimension(:,:), pointer :: cellsOnEdge
+
+ real (kind=RKIND) :: coef
+ real (kind=RKIND), dimension(:), pointer :: dcEdge, dvEdge, areaCell
+ real (kind=RKIND), dimension(:,:), allocatable :: drhoTopOfCell, du2TopOfCell, &
+ drhoTopOfEdge, du2TopOfEdge
+
+ err = 0
+
+ if(.not.richViscOn .and. .not.richDiffOn) return
+
+ nVertLevels = grid % nVertLevels
+ nCells = grid % nCells
+ nEdges = grid % nEdges
+
+ maxLevelCell => grid % maxLevelCell % array
+ maxLevelEdgeTop => grid % maxLevelEdgeTop % array
+ maxLevelEdgeBot => grid % maxLevelEdgeBot % array
+ cellsOnEdge => grid % cellsOnEdge % array
+ dvEdge => grid % dvEdge % array
+ dcEdge => grid % dcEdge % array
+ areaCell => grid % areaCell % array
+
+ allocate( &
+ drhoTopOfCell(nVertLevels+1,nCells+1), drhoTopOfEdge(nVertLevels+1,nEdges+1), &
+ du2TopOfCell(nVertLevels+1,nCells+1), du2TopOfEdge(nVertLevels+1,nEdges+1))
+
+ ! compute density of parcel displaced to next deeper z-level,
+ ! in state % rhoDisplaced
+!maltrud make sure rho is current--check this for redundancy
+! call OcnEquationOfStateRho(grid, 'relative', 0, indexT, indexS, &
+! tracers, rho, err)
+ ! mrp 110324 In order to visualize rhoDisplaced, include the following
+! call OcnEquationOfStateRho(grid, 'relative', 1, indexT, indexS, &
+! tracers, rhoDisplaced, err)
+
+
+ ! drhoTopOfCell(k) = $\rho^*_{k-1}-\rho^*_k$
+ drhoTopOfCell = 0.0
+ do iCell=1,nCells
+ do k=2,maxLevelCell(iCell)
+ drhoTopOfCell(k,iCell) = rho(k-1,iCell) - rhoDisplaced(k-1,iCell)
+ end do
+ end do
+
+ ! interpolate drhoTopOfCell to drhoTopOfEdge
+ do iEdge=1,nEdges
+ cell1 = cellsOnEdge(1,iEdge)
+ cell2 = cellsOnEdge(2,iEdge)
+ do k=2,maxLevelEdgeTop(iEdge)
+ drhoTopOfEdge(k,iEdge) = &
+ (drhoTopOfCell(k,cell1) + &
+ drhoTopOfCell(k,cell2))/2
+ end do
+ end do
+
+ ! du2TopOfEdge(k) = $u_{k-1}-u_k$
+ du2TopOfEdge=0.0
+ do iEdge=1,nEdges
+ do k=2,maxLevelEdgeTop(iEdge)
+ du2TopOfEdge(k,iEdge) = (u(k-1,iEdge) - u(k,iEdge))**2
+ end do
+ end do
+
+ ! interpolate du2TopOfEdge to du2TopOfCell
+ du2TopOfCell(:,:) = 0.0
+ do iEdge=1,nEdges
+ cell1 = cellsOnEdge(1,iEdge)
+ cell2 = cellsOnEdge(2,iEdge)
+ do k=2,maxLevelEdgeBot(iEdge)
+ du2TopOfCell(k,cell1) = du2TopOfCell(k,cell1) &
+ + 0.5 * dcEdge(iEdge) * dvEdge(iEdge) * du2TopOfEdge(k,iEdge)
+ du2TopOfCell(k,cell2) = du2TopOfCell(k,cell2) &
+ + 0.5 * dcEdge(iEdge) * dvEdge(iEdge) * du2TopOfEdge(k,iEdge)
+ end do
+ end do
+ do iCell = 1,nCells
+ do k = 2,maxLevelCell(iCell)
+ du2TopOfCell(k,iCell) = du2TopOfCell(k,iCell) / areaCell(iCell)
+ end do
+ end do
+
+ ! compute RiTopOfEdge using drhoTopOfEdge and du2TopOfEdge
+ ! coef = -g/rho_0/2
+ RiTopOfEdge = 0.0
+ coef = -gravity/1000.0/2.0
+ do iEdge = 1,nEdges
+ do k = 2,maxLevelEdgeTop(iEdge)
+ RiTopOfEdge(k,iEdge) = coef*drhoTopOfEdge(k,iEdge) &
+ *(h_edge(k-1,iEdge)+h_edge(k,iEdge)) &
+ / (du2TopOfEdge(k,iEdge) + 1e-20)
+ end do
+ end do
+
+ ! compute RiTopOfCell using drhoTopOfCell and du2TopOfCell
+ ! coef = -g/rho_0/2
+ RiTopOfCell = 0.0
+ coef = -gravity/1000.0/2.0
+ do iCell = 1,nCells
+ do k = 2,maxLevelCell(iCell)
+ RiTopOfCell(k,iCell) = coef*drhoTopOfCell(k,iCell) &
+ *(h(k-1,iCell)+h(k,iCell)) &
+ / (du2TopOfCell(k,iCell) + 1e-20)
+ end do
+ end do
+
+ deallocate(drhoTopOfCell, drhoTopOfEdge, &
+ du2TopOfCell, du2TopOfEdge)
+
+ !--------------------------------------------------------------------
+
+ end subroutine OcnVmixGetRichNumbers!}}}
+
+!***********************************************************************
+!
+! routine OcnVmixCoefsRichInit
+!
+!> \brief Initializes ocean momentum vertical mixing quantities
+!> \author Doug Jacobsen
+!> \date 19 September 2011
+!> \version SVN:$Id$
+!> \details
+!> This routine initializes a variety of quantities related to
+!> vertical velocity mixing in the ocean. Since a variety of
+!> parameterizations are available, this routine primarily calls the
+!> individual init routines for each parameterization.
+!
+!-----------------------------------------------------------------------
+
+
+ subroutine OcnVmixCoefsRichInit(err)!{{{
+
+ !--------------------------------------------------------------------
+
+ !-----------------------------------------------------------------
+ !
+ ! call individual init routines for each parameterization
+ !
+ !-----------------------------------------------------------------
+
+ integer, intent(out) :: err
+
+ err = 0
+
+ richViscOn = .false.
+ richDiffOn = .false.
+
+ if (config_vert_visc_type.eq.'rich') then
+ richViscOn = .true.
+ endif
+
+ if (config_vert_diff_type.eq.'rich') then
+ richDiffOn = .true.
+ endif
+
+
+ !--------------------------------------------------------------------
+
+ end subroutine OcnVmixCoefsRichInit!}}}
+
+ !!--- TEMPORARY EQUATION OF STATE ROUTINES
+ subroutine equation_of_state(s, grid, k_displaced, displacement_type)!{{{
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ ! This module contains routines necessary for computing the density
+ ! from model temperature and salinity using an equation of state.
+ !
+ ! Input: grid - grid metadata
+ ! s - state: tracers
+ ! k_displaced
+ ! If k_displaced<=0, state % rho is returned with no displaced
+ ! If k_displaced>0,the state % rhoDisplaced is returned, and is for
+ ! a parcel adiabatically displaced from its original level to level
+ ! k_displaced. This does not effect the linear EOS.
+ !
+ ! Output: s - state: computed density
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ implicit none
+
+ type (state_type), intent(inout) :: s
+ type (mesh_type), intent(in) :: grid
+ integer :: k_displaced
+ character(len=8), intent(in) :: displacement_type
+
+ integer, dimension(:), pointer :: maxLevelCell
+ real (kind=RKIND), dimension(:,:), pointer :: rho
+ real (kind=RKIND), dimension(:,:,:), pointer :: tracers
+ integer :: nCells, iCell, k
+ type (dm_info) :: dminfo
+
+ call timer_start("equation_of_state")
+
+ if (config_eos_type.eq.'linear') then
+
+ rho => s % rho % array
+ tracers => s % tracers % array
+ maxLevelCell => grid % maxLevelCell % array
+ nCells = grid % nCells
+
+ do iCell=1,nCells
+ do k=1,maxLevelCell(iCell)
+ ! Linear equation of state
+ rho(k,iCell) = 1000.0*( 1.0 &
+ - 2.5e-4*tracers(s % index_temperature,k,iCell) &
+ + 7.6e-4*tracers(s % index_salinity,k,iCell))
+ end do
+ end do
+
+ elseif (config_eos_type.eq.'jm') then
+
+ call equation_of_state_jm(s, grid, k_displaced, displacement_type)
+
+ else
+ print *, ' Incorrect choice of config_eos_type:',&
+ config_eos_type
+ call dmpar_abort(dminfo)
+ endif
+
+ call timer_stop("equation_of_state")
+
+ end subroutine equation_of_state!}}}
+
+ subroutine equation_of_state_jm(s, grid, k_displaced, displacement_type)!{{{
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ ! This module contains routines necessary for computing the density
+ ! from model temperature and salinity using an equation of state.
+ !
+ ! The UNESCO equation of state computed using the
+ ! potential-temperature-based bulk modulus from Jackett and
+ ! McDougall, JTECH, Vol.12, pp 381-389, April, 1995.
+ !
+ ! Input: grid - grid metadata
+ ! s - state: tracers
+ ! k_displaced
+
+ ! If k_displaced<=0, density is returned with no displaced
+ ! If k_displaced>0,the density returned is that for a parcel
+ ! adiabatically displaced from its original level to level
+ ! k_displaced.
+
+ !
+ ! Output: s - state: computed density
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+ implicit none
+
+ type (state_type), intent(in) :: s
+ type (mesh_type), intent(in) :: grid
+ integer :: k_displaced
+ character(len=8), intent(in) :: displacement_type
+
+ type (dm_info) :: dminfo
+ integer :: iEdge, iCell, iVertex, k
+
+ integer :: nCells, nEdges, nVertices, nVertLevels
+
+
+ real (kind=RKIND), dimension(:), pointer :: &
+ zMidZLevel, pRefEOS
+ real (kind=RKIND), dimension(:,:), pointer :: &
+ rhoPointer
+ real (kind=RKIND), dimension(:,:,:), pointer :: tracers
+
+ integer, dimension(:), pointer :: maxLevelCell
+
+ real (kind=RKIND) :: &
+ TQ,SQ, &! adjusted T,S
+ BULK_MOD, &! Bulk modulus
+ RHO_S, &! density at the surface
+ DRDT0, &! d(density)/d(temperature), for surface
+ DRDS0, &! d(density)/d(salinity ), for surface
+ DKDT, &! d(bulk modulus)/d(pot. temp.)
+ DKDS, &! d(bulk modulus)/d(salinity )
+ SQR,DENOMK, &! work arrays
+ WORK1, WORK2, WORK3, WORK4, T2, depth
+
+ real (kind=RKIND) :: &
+ tmin, tmax, &! valid temperature range for level k
+ smin, smax ! valid salinity range for level k
+
+ real (kind=RKIND), dimension(:), allocatable :: &
+ p, p2 ! temporary pressure scalars
+
+!-----------------------------------------------------------------------
+!
+! UNESCO EOS constants and JMcD bulk modulus constants
+!
+!-----------------------------------------------------------------------
+
+ !*** for density of fresh water (standard UNESCO)
+
+ real (kind=RKIND), parameter :: &
+ unt0 = 999.842594, &
+ unt1 = 6.793952e-2, &
+ unt2 = -9.095290e-3, &
+ unt3 = 1.001685e-4, &
+ unt4 = -1.120083e-6, &
+ unt5 = 6.536332e-9
+
+ !*** for dependence of surface density on salinity (UNESCO)
+
+ real (kind=RKIND), parameter :: &
+ uns1t0 = 0.824493 , &
+ uns1t1 = -4.0899e-3, &
+ uns1t2 = 7.6438e-5, &
+ uns1t3 = -8.2467e-7, &
+ uns1t4 = 5.3875e-9, &
+ unsqt0 = -5.72466e-3, &
+ unsqt1 = 1.0227e-4, &
+ unsqt2 = -1.6546e-6, &
+ uns2t0 = 4.8314e-4
+
+ !*** from Table A1 of Jackett and McDougall
+
+ real (kind=RKIND), parameter :: &
+ bup0s0t0 = 1.965933e+4, &
+ bup0s0t1 = 1.444304e+2, &
+ bup0s0t2 = -1.706103 , &
+ bup0s0t3 = 9.648704e-3, &
+ bup0s0t4 = -4.190253e-5
+
+ real (kind=RKIND), parameter :: &
+ bup0s1t0 = 5.284855e+1, &
+ bup0s1t1 = -3.101089e-1, &
+ bup0s1t2 = 6.283263e-3, &
+ bup0s1t3 = -5.084188e-5
+
+ real (kind=RKIND), parameter :: &
+ bup0sqt0 = 3.886640e-1, &
+ bup0sqt1 = 9.085835e-3, &
+ bup0sqt2 = -4.619924e-4
+
+ real (kind=RKIND), parameter :: &
+ bup1s0t0 = 3.186519 , &
+ bup1s0t1 = 2.212276e-2, &
+ bup1s0t2 = -2.984642e-4, &
+ bup1s0t3 = 1.956415e-6
+
+ real (kind=RKIND), parameter :: &
+ bup1s1t0 = 6.704388e-3, &
+ bup1s1t1 = -1.847318e-4, &
+ bup1s1t2 = 2.059331e-7, &
+ bup1sqt0 = 1.480266e-4
+
+ real (kind=RKIND), parameter :: &
+ bup2s0t0 = 2.102898e-4, &
+ bup2s0t1 = -1.202016e-5, &
+ bup2s0t2 = 1.394680e-7, &
+ bup2s1t0 = -2.040237e-6, &
+ bup2s1t1 = 6.128773e-8, &
+ bup2s1t2 = 6.207323e-10
+
+ integer :: k_test, k_ref
+
+ call timer_start("equation_of_state_jm")
+
+ tracers => s % tracers % array
+
+ nCells = grid % nCells
+ maxLevelCell => grid % maxLevelCell % array
+ nVertLevels = grid % nVertLevels
+ zMidZLevel => grid % zMidZLevel % array
+
+
+! Jackett and McDougall
+ tmin = -2.0 ! valid pot. temp. range
+ tmax = 40.0
+ smin = 0.0 ! valid salinity, in psu
+ smax = 42.0
+
+ ! This could be put in a startup routine.
+ ! Note I am using zMidZlevel, so pressure on top level does
+ ! not include SSH contribution. I am not sure if that matters.
+
+! This function computes pressure in bars from depth in meters
+! using a mean density derived from depth-dependent global
+! average temperatures and salinities from Levitus 1994, and
+! integrating using hydrostatic balance.
+
+ allocate(pRefEOS(nVertLevels),p(nVertLevels),p2(nVertLevels))
+ do k = 1,nVertLevels
+ depth = -zMidZLevel(k)
+ pRefEOS(k) = 0.059808*(exp(-0.025*depth) - 1.0) &
+ + 0.100766*depth + 2.28405e-7*depth**2
+ enddo
+
+ ! If k_displaced=0, in-situ density is returned (no displacement)
+ ! If k_displaced/=0, potential density is returned
+
+ ! if displacement_type = 'relative', potential density is calculated
+ ! referenced to level k + k_displaced
+ ! if displacement_type = 'absolute', potential density is calculated
+ ! referenced to level k_displaced for all k
+ ! NOTE: k_displaced = 0 or > nVertLevels is incompatible with 'absolute'
+ ! so abort if necessary
+
+ if (displacement_type == 'absolute' .and. &
+ (k_displaced <= 0 .or. k_displaced > nVertLevels) ) then
+ write(0,*) 'Abort: In equation_of_state_jm', &
+ ' k_displaced must be between 1 and nVertLevels for ', &
+ 'displacement_type = absolute'
+ call dmpar_abort(dminfo)
+ endif
+
+ if (k_displaced == 0) then
+ rhoPointer => s % rho % array
+ do k=1,nVertLevels
+ p(k) = pRefEOS(k)
+ p2(k) = p(k)*p(k)
+ enddo
+ else ! k_displaced /= 0
+ rhoPointer => s % rhoDisplaced % array
+ do k=1,nVertLevels
+ if (displacement_type == 'relative') then
+ k_test = min(k + k_displaced, nVertLevels)
+ k_ref = max(k_test, 1)
+ else
+ k_test = min(k_displaced, nVertLevels)
+ k_ref = max(k_test, 1)
+ endif
+ p(k) = pRefEOS(k_ref)
+ p2(k) = p(k)*p(k)
+ enddo
+ endif
+
+ do iCell=1,nCells
+ do k=1,maxLevelCell(iCell)
+
+ SQ = max(min(tracers(s%index_salinity,k,iCell),smax),smin)
+ TQ = max(min(tracers(s%index_temperature,k,iCell),tmax),tmin)
+
+ SQR = sqrt(SQ)
+ T2 = TQ*TQ
+
+ !***
+ !*** first calculate surface (p=0) values from UNESCO eqns.
+ !***
+
+ WORK1 = uns1t0 + uns1t1*TQ + &
+ (uns1t2 + uns1t3*TQ + uns1t4*T2)*T2
+ WORK2 = SQR*(unsqt0 + unsqt1*TQ + unsqt2*T2)
+
+ RHO_S = unt1*TQ + (unt2 + unt3*TQ + (unt4 + unt5*TQ)*T2)*T2 &
+ + (uns2t0*SQ + WORK1 + WORK2)*SQ
+
+ !***
+ !*** now calculate bulk modulus at pressure p from
+ !*** Jackett and McDougall formula
+ !***
+
+ WORK3 = bup0s1t0 + bup0s1t1*TQ + &
+ (bup0s1t2 + bup0s1t3*TQ)*T2 + &
+ p(k) *(bup1s1t0 + bup1s1t1*TQ + bup1s1t2*T2) + &
+ p2(k)*(bup2s1t0 + bup2s1t1*TQ + bup2s1t2*T2)
+ WORK4 = SQR*(bup0sqt0 + bup0sqt1*TQ + bup0sqt2*T2 + &
+ bup1sqt0*p(k))
+
+ BULK_MOD = bup0s0t0 + bup0s0t1*TQ + &
+ (bup0s0t2 + bup0s0t3*TQ + bup0s0t4*T2)*T2 + &
+ p(k) *(bup1s0t0 + bup1s0t1*TQ + &
+ (bup1s0t2 + bup1s0t3*TQ)*T2) + &
+ p2(k)*(bup2s0t0 + bup2s0t1*TQ + bup2s0t2*T2) + &
+ SQ*(WORK3 + WORK4)
+
+ DENOMK = 1.0/(BULK_MOD - p(k))
+
+ rhoPointer(k,iCell) = (unt0 + RHO_S)*BULK_MOD*DENOMK
+
+ end do
+ end do
+
+ deallocate(pRefEOS,p,p2)
+
+ call timer_stop("equation_of_state_jm")
+
+ end subroutine equation_of_state_jm!}}}
+
+!***********************************************************************
+
+end module OcnVmixCoefsRich
+
+!|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
+
+!vim: foldmethod=marker
Added: branches/ocean_projects/performance/src/core_ocean/module_OcnVmixCoefsTanh.F
===================================================================
--- branches/ocean_projects/performance/src/core_ocean/module_OcnVmixCoefsTanh.F (rev 0)
+++ branches/ocean_projects/performance/src/core_ocean/module_OcnVmixCoefsTanh.F 2011-09-20 20:17:42 UTC (rev 1013)
@@ -0,0 +1,332 @@
+!|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
+!
+! OcnVmixCoefsTanh
+!
+!> \brief MPAS ocean vertical mixing coefficients
+!> \author Doug Jacobsen
+!> \date 19 September 2011
+!> \version SVN:$Id:$
+!> \details
+!> This module contains the routines for computing
+!> tanhant vertical mixing coefficients.
+!>
+!
+!-----------------------------------------------------------------------
+
+module OcnVmixCoefsTanh
+
+ use grid_types
+ use configure
+ use timer
+
+ implicit none
+ private
+ save
+
+ !--------------------------------------------------------------------
+ !
+ ! Public parameters
+ !
+ !--------------------------------------------------------------------
+
+ !--------------------------------------------------------------------
+ !
+ ! Public member functions
+ !
+ !--------------------------------------------------------------------
+
+ private :: OcnVelVmixCoefsTanh, &
+ OcnTracerVmixCoefsTanh
+
+ public :: OcnVmixCoefsTanhBuild, &
+ OcnVmixCoefsTanhInit
+
+ !--------------------------------------------------------------------
+ !
+ ! Private module variables
+ !
+ !--------------------------------------------------------------------
+
+ logical :: tanhViscOn, tanhDiffOn
+
+!***********************************************************************
+
+contains
+
+!***********************************************************************
+!
+! routine OcnVmixCoefsTanhBuild
+!
+!> \brief Computes coefficients for vertical mixing
+!> \author Doug Jacobsen
+!> \date 19 September 2011
+!> \version SVN:$Id$
+!> \details
+!> This routine computes the vertical mixing coefficients for momentum
+!> and tracers based user choices of mixing parameterization.
+!
+!-----------------------------------------------------------------------
+
+ subroutine OcnVmixCoefsTanhBuild(grid, s, d, err)!{{{
+
+ !-----------------------------------------------------------------
+ !
+ ! input variables
+ !
+ !-----------------------------------------------------------------
+
+ type (mesh_type), intent(in) :: &
+ grid !< Input: grid information
+
+ !-----------------------------------------------------------------
+ !
+ ! input/output variables
+ !
+ !-----------------------------------------------------------------
+
+ type (state_type), intent(inout) :: &
+ s !< Input/Output: state information
+
+ type (diagnostics_type), intent(inout) :: &
+ d !< Input/Output: diagnostic information
+
+ !-----------------------------------------------------------------
+ !
+ ! output variables
+ !
+ !-----------------------------------------------------------------
+
+ integer, intent(out) :: err
+
+ !-----------------------------------------------------------------
+ !
+ ! local variables
+ !
+ !-----------------------------------------------------------------
+
+ integer :: err1, err2
+
+ real (kind=RKIND), dimension(:,:), pointer :: &
+ vertViscTopOfEdge, vertDiffTopOfCell
+
+ !-----------------------------------------------------------------
+ !
+ ! call relevant routines for computing tendencies
+ ! note that the user can choose multiple options and the
+ ! tendencies will be added together
+ !
+ !-----------------------------------------------------------------
+
+ err = 0
+ if((.not.tanhViscOn) .and. (.not.tanhDiffOn)) return
+
+ vertViscTopOfEdge => d % vertViscTopOfEdge % array
+ vertDiffTopOfCell => d % vertDiffTopOfCell % array
+
+ call OcnVelVmixCoefsTanh(grid, vertViscTopOfEdge, err1)
+ call OcnTracerVmixCoefsTanh(grid, vertDiffTopOfCell, err2)
+
+ err = err1 .or. err2
+
+ !--------------------------------------------------------------------
+
+ end subroutine OcnVmixCoefsTanhBuild!}}}
+
+!***********************************************************************
+!
+! routine OcnVelVmixCoefsTanh
+!
+!> \brief Computes coefficients for vertical momentum mixing
+!> \author Doug Jacobsen
+!> \date 19 September 2011
+!> \version SVN:$Id$
+!> \details
+!> This routine computes the tanh vertical mixing coefficients for momentum
+!
+!-----------------------------------------------------------------------
+
+ subroutine OcnVelVmixCoefsTanh(grid, vertViscTopOfEdge, err)!{{{
+
+ !-----------------------------------------------------------------
+ !
+ ! input variables
+ !
+ !-----------------------------------------------------------------
+
+ type (mesh_type), intent(in) :: &
+ grid !< Input: grid information
+
+ !-----------------------------------------------------------------
+ !
+ ! input/output variables
+ !
+ !-----------------------------------------------------------------
+
+ !-----------------------------------------------------------------
+ !
+ ! output variables
+ !
+ !-----------------------------------------------------------------
+
+ real (kind=RKIND), dimension(:,:), intent(out) :: vertViscTopOfEdge
+
+ integer, intent(out) :: err
+
+ !-----------------------------------------------------------------
+ !
+ ! local variables
+ !
+ !-----------------------------------------------------------------
+
+ integer :: k, nVertLevels
+
+ real (kind=RKIND), dimension(:), pointer :: zTopZLevel
+
+ err = 0
+
+ if(.not.tanhViscOn) return
+
+ nVertLevels = grid % nVertLevels
+ zTopZLevel => grid % zTopZLevel % array
+
+ do k=1,nVertLevels+1
+ vertViscTopOfEdge(k,:) = -(config_max_visc_tanh-config_min_visc_tanh)/2.0 &
+ *tanh(-(zTopZLevel(k)-config_ZMid_tanh) &
+ /config_zWidth_tanh) &
+ + (config_max_visc_tanh+config_min_visc_tanh)/2
+ end do
+
+
+ !--------------------------------------------------------------------
+
+ end subroutine OcnVelVmixCoefsTanh!}}}
+
+!***********************************************************************
+!
+! routine OcnTracerVmixCoefsTanh
+!
+!> \brief Computes coefficients for vertical tracer mixing
+!> \author Doug Jacobsen
+!> \date 19 September 2011
+!> \version SVN:$Id$
+!> \details
+!> This routine computes the tanh vertical mixing coefficients for tracers
+!
+!-----------------------------------------------------------------------
+
+ subroutine OcnTracerVmixCoefsTanh(grid, vertDiffTopOfCell, err)!{{{
+
+ !-----------------------------------------------------------------
+ !
+ ! input variables
+ !
+ !-----------------------------------------------------------------
+
+ type (mesh_type), intent(in) :: &
+ grid !< Input: grid information
+
+ !-----------------------------------------------------------------
+ !
+ ! input/output variables
+ !
+ !-----------------------------------------------------------------
+
+ !-----------------------------------------------------------------
+ !
+ ! output variables
+ !
+ !-----------------------------------------------------------------
+
+ real (kind=RKIND), dimension(:,:), intent(out) :: vertDiffTopOfCell
+
+ integer, intent(out) :: err
+
+ !-----------------------------------------------------------------
+ !
+ ! local variables
+ !
+ !-----------------------------------------------------------------
+
+ integer :: k, nVertLevels
+
+ real (kind=RKIND), dimension(:), pointer :: zTopZLevel
+
+ err = 0
+
+ if(.not.tanhDiffOn) return
+
+ nVertLevels = grid % nVertLevels
+ zTopZLevel => grid % zTopZLevel % array
+
+ do k=1,nVertLevels+1
+ vertDiffTopOfCell(k,:) = -(config_max_diff_tanh-config_min_diff_tanh)/2.0 &
+ *tanh(-(zTopZLevel(k)-config_ZMid_tanh) &
+ /config_zWidth_tanh) &
+ + (config_max_diff_tanh+config_min_diff_tanh)/2
+ end do
+
+
+ !--------------------------------------------------------------------
+
+ end subroutine OcnTracerVmixCoefsTanh!}}}
+
+
+!***********************************************************************
+!
+! routine OcnVmixCoefsTanhInit
+!
+!> \brief Initializes ocean vertical mixing quantities
+!> \author Doug Jacobsen
+!> \date 19 September 2011
+!> \version SVN:$Id$
+!> \details
+!> This routine initializes a variety of quantities related to
+!> tanh vertical mixing in the ocean.
+!
+!-----------------------------------------------------------------------
+
+
+ subroutine OcnVmixCoefsTanhInit(err)!{{{
+
+ !--------------------------------------------------------------------
+
+ !-----------------------------------------------------------------
+ !
+ ! call individual init routines for each parameterization
+ !
+ !-----------------------------------------------------------------
+
+ integer, intent(out) :: err
+
+ err = 0
+
+ tanhViscOn = .false.
+ tanhDiffOn = .false.
+
+ if (config_vert_visc_type.eq.'tanh') then
+ tanhViscOn = .true.
+ endif
+
+ if (config_vert_diff_type.eq.'tanh') then
+ tanhDiffOn = .true.
+ endif
+
+ if(tanhViscOn .or. tanhDiffOn) then
+ if (config_vert_grid_type.ne.'zlevel') then
+ write(0,*) 'Abort: config_vert_diff_type.eq.tanh may only', &
+ ' use config_vert_grid_type of zlevel at this time'
+ err = 1
+ endif
+ endif
+
+ !--------------------------------------------------------------------
+
+ end subroutine OcnVmixCoefsTanhInit!}}}
+
+!***********************************************************************
+
+end module OcnVmixCoefsTanh
+
+!|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
+
+!vim: foldmethod=marker
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-20 17:15:08 UTC (rev 1012)
+++ branches/ocean_projects/performance/src/core_ocean/module_mpas_core.F 2011-09-20 20:17:42 UTC (rev 1013)
@@ -16,6 +16,7 @@
use OcnRestoring
use OcnEquationOfState
+ use OcnVmix
type (io_output_object) :: restart_obj
integer :: restart_frame
@@ -97,7 +98,12 @@
call OcnRestoringInit(err)
call OcnEquationOfStateInit(err)
+ call OcnVmixInit(err)
+ if(err) then
+ call dmpar_abort(dminfo)
+ endif
+
! 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-20 17:15:08 UTC (rev 1012)
+++ branches/ocean_projects/performance/src/core_ocean/module_time_integration.F 2011-09-20 20:17:42 UTC (rev 1013)
@@ -23,6 +23,7 @@
use OcnRestoring
use OcnEquationOfState
+ use OcnVmix
contains
@@ -88,7 +89,7 @@
type (domain_type), intent(inout) :: domain
real (kind=RKIND), intent(in) :: dt
- integer :: iCell, k, i
+ integer :: iCell, k, i, err
type (block_type), pointer :: block
type (state_type) :: provis
@@ -179,7 +180,7 @@
block => domain % blocklist
do while (associated(block))
if (.not.config_implicit_vertical_mix) then
- call compute_vertical_mix_coefficients(provis, block % diagnostics, block % mesh)
+ call OcnVmixCoefs(block%mesh, provis, block%diagnostics, err)
end if
call compute_tend_h(block % tend, provis, block % diagnostics, block % mesh)
call compute_tend_u(block % tend, provis, block % diagnostics, block % mesh)
@@ -308,49 +309,12 @@
if (config_implicit_vertical_mix) then
call timer_start("RK4-implicit vert mix")
- allocate(A(nVertLevels),C(nVertLevels),uTemp(nVertLevels), &
- tracersTemp(num_tracers,nVertLevels))
- call compute_vertical_mix_coefficients(block % state % time_levs(2) % state, block % diagnostics, block % mesh)
+ call OcnVmixCoefs(block%mesh, provis, block%diagnostics, err)
- !
- ! Implicit vertical solve for momentum
- !
- do iEdge=1,nEdges
- if (maxLevelEdgeTop(iEdge).gt.0) then
-
- ! Compute A(k), C(k) for momentum
- ! mrp 110315 efficiency note: for z-level, could precompute
- ! -2.0*dt/(h(k)_h(k+1))/h(k) in setup
- ! h_edge is computed in compute_solve_diag, and is not available yet.
- ! This could be removed if hZLevel used instead.
- cell1 = block % mesh % cellsOnEdge % array(1,iEdge)
- cell2 = block % mesh % cellsOnEdge % array(2,iEdge)
- do k=1,maxLevelEdgeTop(iEdge)
- h_edge(k,iEdge) = 0.5 * (h(k,cell1) + h(k,cell2))
- end do
-
- do k=1,maxLevelEdgeTop(iEdge)-1
- A(k) = -2.0*dt*vertViscTopOfEdge(k+1,iEdge) &
- / (h_edge(k,iEdge) + h_edge(k+1,iEdge)) &
- / h_edge(k,iEdge)
- enddo
- A(maxLevelEdgeTop(iEdge)) = -dt*config_bottom_drag_coeff &
- *sqrt(2.0*ke_edge(k,iEdge))/h_edge(k,iEdge)
-
- C(1) = 1 - A(1)
- do k=2,maxLevelEdgeTop(iEdge)
- C(k) = 1 - A(k) - A(k-1)
- enddo
-
- call tridiagonal_solve(A,C,A,u(:,iEdge),uTemp,maxLevelEdgeTop(iEdge))
-
- u(1:maxLevelEdgeTop(iEdge),iEdge) = uTemp(1:maxLevelEdgeTop(iEdge))
- u(maxLevelEdgeTop(iEdge)+1:nVertLevels,iEdge) = 0.0
-
- end if
- end do
-
+ call OcnVelVmixTendImplicit(block%mesh, dt, ke_edge, vertViscTopOfEdge,h, &
+ h_edge, u, err)
+
! mrp 110718 filter btr mode out of u
if (config_rk_filter_btr_mode) then
call filter_btr_mode_u(block % state % time_levs(2) % state, block % mesh)
@@ -360,32 +324,9 @@
!
! Implicit vertical solve for tracers
!
- do iCell=1,nCells
- ! Compute A(k), C(k) for tracers
- ! mrp 110315 efficiency note: for z-level, could precompute
- ! -2.0*dt/(h(k)_h(k+1))/h(k) in setup
- do k=1,maxLevelCell(iCell)-1
- A(k) = -2.0*dt*vertDiffTopOfCell(k+1,iCell) &
- / (h(k,iCell) + h(k+1,iCell)) &
- / h(k,iCell)
- enddo
-
- A(maxLevelCell(iCell)) = 0.0
-
- C(1) = 1 - A(1)
- do k=2,maxLevelCell(iCell)
- C(k) = 1 - A(k) - A(k-1)
- enddo
-
- call tridiagonal_solve_mult(A,C,A,tracers(:,:,iCell), &
- tracersTemp, maxLevelCell(iCell), &
- nVertLevels,num_tracers)
-
- tracers(:,1:maxLevelCell(iCell),iCell) = tracersTemp(:,1:maxLevelCell(iCell))
- tracers(:,maxLevelCell(iCell)+1:nVertLevels,iCell) = -1e34
- end do
- deallocate(A,C,uTemp,tracersTemp)
- call timer_stop("RK4-implicit vert mix")
+
+ call OcnTracerVmixTendImplicit(block%mesh, dt, vertDiffTopOfCell, h, &
+ tracers, err)
end if
! mrp 110725 momentum decay term
@@ -1857,29 +1798,8 @@
!
! velocity tendency: vertical mixing d/dz( nu_v du/dz))
!
- if (.not.config_implicit_vertical_mix) then
- call timer_start("compute_tend_u-explicit vert mix")
- allocate(fluxVertTop(nVertLevels+1))
- fluxVertTop(1) = 0.0
- do iEdge=1,grid % nEdgesSolve
-
- do k=2,maxLevelEdgeTop(iEdge)
- fluxVertTop(k) = vertViscTopOfEdge(k,iEdge) &
- * ( u(k-1,iEdge) - u(k,iEdge) ) &
- * 2 / (h_edge(k-1,iEdge) + h_edge(k,iEdge))
- enddo
- fluxVertTop(maxLevelEdgeTop(iEdge)+1) = 0.0
+ call OcnVelVmixTendExplicit(grid, u, h_edge, vertViscTopOfEdge, tend_u, err)
- do k=1,maxLevelEdgeTop(iEdge)
- tend_u(k,iEdge) = tend_u(k,iEdge) &
- + (fluxVertTop(k) - fluxVertTop(k+1)) &
- / h_edge(k,iEdge)
- enddo
-
- end do
- deallocate(fluxVertTop)
- call timer_stop("compute_tend_u-explicit vert mix")
- endif
call timer_stop("compute_tend_u")
end subroutine compute_tend_u!}}}
@@ -2379,40 +2299,9 @@
!
! tracer tendency: vertical diffusion h d/dz( \kappa_v d\phi/dz)
!
- if (.not.config_implicit_vertical_mix) then
- call timer_start("compute_scalar_tend-explicit vert diff")
- allocate(fluxVertTop(num_tracers,nVertLevels+1))
- fluxVertTop(:,1) = 0.0
- do iCell=1,nCellsSolve
+ call OcnTracerVmixTendExplicit(grid, h, vertDiffTopOfCell, tracers, tend_tr, err)
- do k=2,maxLevelCell(iCell)
- do iTracer=1,num_tracers
- ! compute \kappa_v d\phi/dz
- fluxVertTop(iTracer,k) = vertDiffTopOfCell(k,iCell) &
- * ( tracers(iTracer,k-1,iCell) &
- - tracers(iTracer,k ,iCell) ) &
- * 2 / (h(k-1,iCell) + h(k,iCell))
- enddo
- enddo
- fluxVertTop(:,maxLevelCell(iCell)+1) = 0.0
-
- do k=1,maxLevelCell(iCell)
- do iTracer=1,num_tracers
- ! This is h d/dz( fluxVertTop) but h and dz cancel, so
- ! reduces to delta( fluxVertTop)
- tend_tr(iTracer,k,iCell) = tend_tr(iTracer,k,iCell) &
- + fluxVertTop(iTracer,k) - fluxVertTop(iTracer,k+1)
- enddo
- enddo
-!print '(a,50e12.2)', 'fluxVertTop',fluxVertTop(3,1:maxLevelCell(iCell)+1)
-!print '(a,50e12.2)', 'tend_tr ',tend_tr(3,1,1:maxLevelCell(iCell))
-
- enddo ! iCell loop
- deallocate(fluxVertTop)
- call timer_stop("compute_scalar_tend-explicit vert diff")
- endif
-
! mrp 110516 printing
!print *, 'tend_tr 2',minval(tend_tr(3,1,1:nCells)),&
! maxval(tend_tr(3,1,1:nCells))
@@ -2847,12 +2736,18 @@
!
! equation of state
!
- call OcnEquationOfStateRho(grid, 'relative', 0, s%index_temperature, s%index_salinity, &
- tracers, rho, err)
+! call OcnEquationOfStateRho(grid, 'relative', 0, s%index_temperature, s%index_salinity, &
+! tracers, rho, err)
! mrp 110324 In order to visualize rhoDisplaced, include the following
- call OcnEquationOfStateRho(grid, 'relative', 1, s%index_temperature, s%index_salinity, &
- tracers, rhoDisplaced, err)
+! call OcnEquationOfStateRho(grid, 'relative', 1, s%index_temperature, s%index_salinity, &
+! tracers, rhoDisplaced, err)
+ if (config_vert_grid_type.eq.'zlevel') then
+ call equation_of_state(s,grid,0,'relative')
+ ! mrp 110324 In order to visualize rhoDisplaced, include the following
+ call equation_of_state(s, grid, 1,'relative')
+ endif
+
!
! Pressure
! This section must be after computing rho
@@ -3091,12 +2986,15 @@
! compute density of parcel displaced to next deeper z-level,
! in state % rhoDisplaced
!maltrud make sure rho is current--check this for redundancy
- call OcnEquationOfStateRho(grid, 'relative', 0, s%index_temperature, s%index_salinity, &
- tracers, rho, err)
+! call OcnEquationOfStateRho(grid, 'relative', 0, s%index_temperature, s%index_salinity, &
+! tracers, rho, err)
! mrp 110324 In order to visualize rhoDisplaced, include the following
- call OcnEquationOfStateRho(grid, 'relative', 1, s%index_temperature, s%index_salinity, &
- tracers, rhoDisplaced, err)
+! call OcnEquationOfStateRho(grid, 'relative', 1, s%index_temperature, s%index_salinity, &
+! tracers, rhoDisplaced, err)
+ call equation_of_state(s, grid, 0, 'relative')
+ call equation_of_state(s, grid, 1, 'relative')
+
! drhoTopOfCell(k) = $\rho^*_{k-1}-\rho^*_k$
drhoTopOfCell = 0.0
do iCell=1,nCells
@@ -3353,6 +3251,321 @@
end subroutine enforce_boundaryEdge!}}}
+ subroutine equation_of_state(s, grid, k_displaced, displacement_type)!{{{
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ ! This module contains routines necessary for computing the density
+ ! from model temperature and salinity using an equation of state.
+ !
+ ! Input: grid - grid metadata
+ ! s - state: tracers
+ ! k_displaced
+ ! If k_displaced<=0, state % rho is returned with no displaced
+ ! If k_displaced>0,the state % rhoDisplaced is returned, and is for
+ ! a parcel adiabatically displaced from its original level to level
+ ! k_displaced. This does not effect the linear EOS.
+ !
+ ! Output: s - state: computed density
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ implicit none
+
+ type (state_type), intent(inout) :: s
+ type (mesh_type), intent(in) :: grid
+ integer :: k_displaced
+ character(len=8), intent(in) :: displacement_type
+
+ integer, dimension(:), pointer :: maxLevelCell
+ real (kind=RKIND), dimension(:,:), pointer :: rho
+ real (kind=RKIND), dimension(:,:,:), pointer :: tracers
+ integer :: nCells, iCell, k
+ type (dm_info) :: dminfo
+
+ call timer_start("equation_of_state")
+
+ if (config_eos_type.eq.'linear') then
+
+ rho => s % rho % array
+ tracers => s % tracers % array
+ maxLevelCell => grid % maxLevelCell % array
+ nCells = grid % nCells
+
+ do iCell=1,nCells
+ do k=1,maxLevelCell(iCell)
+ ! Linear equation of state
+ rho(k,iCell) = 1000.0*( 1.0 &
+ - 2.5e-4*tracers(s % index_temperature,k,iCell) &
+ + 7.6e-4*tracers(s % index_salinity,k,iCell))
+ end do
+ end do
+
+ elseif (config_eos_type.eq.'jm') then
+
+ call equation_of_state_jm(s, grid, k_displaced, displacement_type)
+
+ else
+ print *, ' Incorrect choice of config_eos_type:',&
+ config_eos_type
+ call dmpar_abort(dminfo)
+ endif
+
+ call timer_stop("equation_of_state")
+
+ end subroutine equation_of_state!}}}
+
+ subroutine equation_of_state_jm(s, grid, k_displaced, displacement_type)!{{{
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ ! This module contains routines necessary for computing the density
+ ! from model temperature and salinity using an equation of state.
+ !
+ ! The UNESCO equation of state computed using the
+ ! potential-temperature-based bulk modulus from Jackett and
+ ! McDougall, JTECH, Vol.12, pp 381-389, April, 1995.
+ !
+ ! Input: grid - grid metadata
+ ! s - state: tracers
+ ! k_displaced
+
+ ! If k_displaced<=0, density is returned with no displaced
+ ! If k_displaced>0,the density returned is that for a parcel
+ ! adiabatically displaced from its original level to level
+ ! k_displaced.
+
+ !
+ ! Output: s - state: computed density
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+ implicit none
+
+ type (state_type), intent(in) :: s
+ type (mesh_type), intent(in) :: grid
+ integer :: k_displaced
+ character(len=8), intent(in) :: displacement_type
+
+ type (dm_info) :: dminfo
+ integer :: iEdge, iCell, iVertex, k
+
+ integer :: nCells, nEdges, nVertices, nVertLevels
+
+
+ real (kind=RKIND), dimension(:), pointer :: &
+ zMidZLevel, pRefEOS
+ real (kind=RKIND), dimension(:,:), pointer :: &
+ rhoPointer
+ real (kind=RKIND), dimension(:,:,:), pointer :: tracers
+
+ integer, dimension(:), pointer :: maxLevelCell
+
+ real (kind=RKIND) :: &
+ TQ,SQ, &! adjusted T,S
+ BULK_MOD, &! Bulk modulus
+ RHO_S, &! density at the surface
+ DRDT0, &! d(density)/d(temperature), for surface
+ DRDS0, &! d(density)/d(salinity ), for surface
+ DKDT, &! d(bulk modulus)/d(pot. temp.)
+ DKDS, &! d(bulk modulus)/d(salinity )
+ SQR,DENOMK, &! work arrays
+ WORK1, WORK2, WORK3, WORK4, T2, depth
+
+ real (kind=RKIND) :: &
+ tmin, tmax, &! valid temperature range for level k
+ smin, smax ! valid salinity range for level k
+
+ real (kind=RKIND), dimension(:), allocatable :: &
+ p, p2 ! temporary pressure scalars
+
+!-----------------------------------------------------------------------
+!
+! UNESCO EOS constants and JMcD bulk modulus constants
+!
+!-----------------------------------------------------------------------
+
+ !*** for density of fresh water (standard UNESCO)
+
+ real (kind=RKIND), parameter :: &
+ unt0 = 999.842594, &
+ unt1 = 6.793952e-2, &
+ unt2 = -9.095290e-3, &
+ unt3 = 1.001685e-4, &
+ unt4 = -1.120083e-6, &
+ unt5 = 6.536332e-9
+
+ !*** for dependence of surface density on salinity (UNESCO)
+
+ real (kind=RKIND), parameter :: &
+ uns1t0 = 0.824493 , &
+ uns1t1 = -4.0899e-3, &
+ uns1t2 = 7.6438e-5, &
+ uns1t3 = -8.2467e-7, &
+ uns1t4 = 5.3875e-9, &
+ unsqt0 = -5.72466e-3, &
+ unsqt1 = 1.0227e-4, &
+ unsqt2 = -1.6546e-6, &
+ uns2t0 = 4.8314e-4
+
+ !*** from Table A1 of Jackett and McDougall
+
+ real (kind=RKIND), parameter :: &
+ bup0s0t0 = 1.965933e+4, &
+ bup0s0t1 = 1.444304e+2, &
+ bup0s0t2 = -1.706103 , &
+ bup0s0t3 = 9.648704e-3, &
+ bup0s0t4 = -4.190253e-5
+
+ real (kind=RKIND), parameter :: &
+ bup0s1t0 = 5.284855e+1, &
+ bup0s1t1 = -3.101089e-1, &
+ bup0s1t2 = 6.283263e-3, &
+ bup0s1t3 = -5.084188e-5
+
+ real (kind=RKIND), parameter :: &
+ bup0sqt0 = 3.886640e-1, &
+ bup0sqt1 = 9.085835e-3, &
+ bup0sqt2 = -4.619924e-4
+
+ real (kind=RKIND), parameter :: &
+ bup1s0t0 = 3.186519 , &
+ bup1s0t1 = 2.212276e-2, &
+ bup1s0t2 = -2.984642e-4, &
+ bup1s0t3 = 1.956415e-6
+
+ real (kind=RKIND), parameter :: &
+ bup1s1t0 = 6.704388e-3, &
+ bup1s1t1 = -1.847318e-4, &
+ bup1s1t2 = 2.059331e-7, &
+ bup1sqt0 = 1.480266e-4
+
+ real (kind=RKIND), parameter :: &
+ bup2s0t0 = 2.102898e-4, &
+ bup2s0t1 = -1.202016e-5, &
+ bup2s0t2 = 1.394680e-7, &
+ bup2s1t0 = -2.040237e-6, &
+ bup2s1t1 = 6.128773e-8, &
+ bup2s1t2 = 6.207323e-10
+
+ integer :: k_test, k_ref
+
+ call timer_start("equation_of_state_jm")
+
+ tracers => s % tracers % array
+
+ nCells = grid % nCells
+ maxLevelCell => grid % maxLevelCell % array
+ nVertLevels = grid % nVertLevels
+ zMidZLevel => grid % zMidZLevel % array
+
+
+! Jackett and McDougall
+ tmin = -2.0 ! valid pot. temp. range
+ tmax = 40.0
+ smin = 0.0 ! valid salinity, in psu
+ smax = 42.0
+
+ ! This could be put in a startup routine.
+ ! Note I am using zMidZlevel, so pressure on top level does
+ ! not include SSH contribution. I am not sure if that matters.
+
+! This function computes pressure in bars from depth in meters
+! using a mean density derived from depth-dependent global
+! average temperatures and salinities from Levitus 1994, and
+! integrating using hydrostatic balance.
+
+ allocate(pRefEOS(nVertLevels),p(nVertLevels),p2(nVertLevels))
+ do k = 1,nVertLevels
+ depth = -zMidZLevel(k)
+ pRefEOS(k) = 0.059808*(exp(-0.025*depth) - 1.0) &
+ + 0.100766*depth + 2.28405e-7*depth**2
+ enddo
+
+ ! If k_displaced=0, in-situ density is returned (no displacement)
+ ! If k_displaced/=0, potential density is returned
+
+ ! if displacement_type = 'relative', potential density is calculated
+ ! referenced to level k + k_displaced
+ ! if displacement_type = 'absolute', potential density is calculated
+ ! referenced to level k_displaced for all k
+ ! NOTE: k_displaced = 0 or > nVertLevels is incompatible with 'absolute'
+ ! so abort if necessary
+
+ if (displacement_type == 'absolute' .and. &
+ (k_displaced <= 0 .or. k_displaced > nVertLevels) ) then
+ write(0,*) 'Abort: In equation_of_state_jm', &
+ ' k_displaced must be between 1 and nVertLevels for ', &
+ 'displacement_type = absolute'
+ call dmpar_abort(dminfo)
+ endif
+
+ if (k_displaced == 0) then
+ rhoPointer => s % rho % array
+ do k=1,nVertLevels
+ p(k) = pRefEOS(k)
+ p2(k) = p(k)*p(k)
+ enddo
+ else ! k_displaced /= 0
+ rhoPointer => s % rhoDisplaced % array
+ do k=1,nVertLevels
+ if (displacement_type == 'relative') then
+ k_test = min(k + k_displaced, nVertLevels)
+ k_ref = max(k_test, 1)
+ else
+ k_test = min(k_displaced, nVertLevels)
+ k_ref = max(k_test, 1)
+ endif
+ p(k) = pRefEOS(k_ref)
+ p2(k) = p(k)*p(k)
+ enddo
+ endif
+
+ do iCell=1,nCells
+ do k=1,maxLevelCell(iCell)
+
+ SQ = max(min(tracers(s%index_salinity,k,iCell),smax),smin)
+ TQ = max(min(tracers(s%index_temperature,k,iCell),tmax),tmin)
+
+ SQR = sqrt(SQ)
+ T2 = TQ*TQ
+
+ !***
+ !*** first calculate surface (p=0) values from UNESCO eqns.
+ !***
+
+ WORK1 = uns1t0 + uns1t1*TQ + &
+ (uns1t2 + uns1t3*TQ + uns1t4*T2)*T2
+ WORK2 = SQR*(unsqt0 + unsqt1*TQ + unsqt2*T2)
+
+ RHO_S = unt1*TQ + (unt2 + unt3*TQ + (unt4 + unt5*TQ)*T2)*T2 &
+ + (uns2t0*SQ + WORK1 + WORK2)*SQ
+
+ !***
+ !*** now calculate bulk modulus at pressure p from
+ !*** Jackett and McDougall formula
+ !***
+
+ WORK3 = bup0s1t0 + bup0s1t1*TQ + &
+ (bup0s1t2 + bup0s1t3*TQ)*T2 + &
+ p(k) *(bup1s1t0 + bup1s1t1*TQ + bup1s1t2*T2) + &
+ p2(k)*(bup2s1t0 + bup2s1t1*TQ + bup2s1t2*T2)
+ WORK4 = SQR*(bup0sqt0 + bup0sqt1*TQ + bup0sqt2*T2 + &
+ bup1sqt0*p(k))
+
+ BULK_MOD = bup0s0t0 + bup0s0t1*TQ + &
+ (bup0s0t2 + bup0s0t3*TQ + bup0s0t4*T2)*T2 + &
+ p(k) *(bup1s0t0 + bup1s0t1*TQ + &
+ (bup1s0t2 + bup1s0t3*TQ)*T2) + &
+ p2(k)*(bup2s0t0 + bup2s0t1*TQ + bup2s0t2*T2) + &
+ SQ*(WORK3 + WORK4)
+
+ DENOMK = 1.0/(BULK_MOD - p(k))
+
+ rhoPointer(k,iCell) = (unt0 + RHO_S)*BULK_MOD*DENOMK
+
+ end do
+ end do
+
+ deallocate(pRefEOS,p,p2)
+
+ call timer_stop("equation_of_state_jm")
+
+ end subroutine equation_of_state_jm!}}}
+
subroutine tridiagonal_solve(a,b,c,r,x,n)!{{{
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Solve the matrix equation Ax=r for x, where A is tridiagonal.
</font>
</pre>