<p><b>dwj07@fsu.edu</b> 2011-09-23 15:26:14 -0600 (Fri, 23 Sep 2011)</p><p><br>
        -- BRANCH COMMIT --<br>
<br>
        Splitting tendency subroutines into their own module to prepare for timestepping modules.<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-23 18:09:34 UTC (rev 1018)
+++ branches/ocean_projects/performance/src/core_ocean/Makefile        2011-09-23 21:26:14 UTC (rev 1019)
@@ -35,6 +35,7 @@
         module_OcnVmixCoefsTanh.o \
         module_OcnRestoring.o \
         module_OcnEquationOfState.o \
+         module_OcnTendency.o \
module_time_integration.o \
module_global_diagnostics.o
@@ -48,39 +49,42 @@
module_advection.o:
-module_time_integration.o: module_OcnThickHadv.o \
-                                                        module_OcnThickVadv.o \
-                                                        module_OcnVelCoriolis.o \
-                                                        module_OcnVelVadv.o \
-                                                        module_OcnVelHmix.o \
-                                                        module_OcnVelHmixDel2.o \
-                                                        module_OcnVelHmixDel4.o \
-                                                        module_OcnVelForcing.o \
-                                                        module_OcnVelForcingWindStress.o \
-                                                        module_OcnVelForcingBottomDrag.o \
-                                                        module_OcnVelPressureGrad.o \
-                                                        module_OcnTracerVadv.o \
-                                                        module_OcnTracerVadvSpline.o \
-                                                        module_OcnTracerVadvSpline2.o \
-                                                        module_OcnTracerVadvSpline3.o \
-                                                        module_OcnTracerVadvStencil.o \
-                                                        module_OcnTracerVadvStencil2.o \
-                                                        module_OcnTracerVadvStencil3.o \
-                                                        module_OcnTracerVadvStencil4.o \
-                                                        module_OcnTracerHadv.o \
-                                                        module_OcnTracerHadv2.o \
-                                                        module_OcnTracerHadv3.o \
-                                                        module_OcnTracerHadv4.o \
-                                                        module_OcnTracerHmix.o \
-                                                        module_OcnTracerHmixDel2.o \
-                                                        module_OcnTracerHmixDel4.o \
-                                                        module_OcnVmix.o \
-                                                        module_OcnVmixCoefsConst.o \
-                                                        module_OcnVmixCoefsRich.o \
-                                                        module_OcnVmixCoefsTanh.o \
-                                                        module_OcnRestoring.o \
-                                                        module_OcnEquationOfState.o
+module_time_integration.o: module_OcnThickHadv.o \
+         module_OcnThickVadv.o \
+                                                 module_OcnVelCoriolis.o \
+                                                 module_OcnVelVadv.o \
+                                                 module_OcnVelHmix.o \
+                                                 module_OcnVelHmixDel2.o \
+                                                 module_OcnVelHmixDel4.o \
+                                                 module_OcnVelForcing.o \
+                                                 module_OcnVelForcingWindStress.o \
+                                                 module_OcnVelForcingBottomDrag.o \
+                                                 module_OcnVelPressureGrad.o \
+                                                 module_OcnTracerVadv.o \
+                                                 module_OcnTracerVadvSpline.o \
+                                                 module_OcnTracerVadvSpline2.o \
+                                                 module_OcnTracerVadvSpline3.o \
+                                                 module_OcnTracerVadvStencil.o \
+                                                 module_OcnTracerVadvStencil2.o \
+                                                 module_OcnTracerVadvStencil3.o \
+                                                 module_OcnTracerVadvStencil4.o \
+                                                 module_OcnTracerHadv.o \
+                                                 module_OcnTracerHadv2.o \
+                                                 module_OcnTracerHadv3.o \
+                                                 module_OcnTracerHadv4.o \
+                                                 module_OcnTracerHmix.o \
+                                                 module_OcnTracerHmixDel2.o \
+                                                 module_OcnTracerHmixDel4.o \
+                                                 module_OcnVmix.o \
+                                                 module_OcnVmixCoefsConst.o \
+                                                 module_OcnVmixCoefsRich.o \
+                                                 module_OcnVmixCoefsTanh.o \
+                                                 module_OcnRestoring.o \
+                                                 module_OcnEquationOfState.o \
+                                                 module_OcnTendency.o
+module_OcnTendency.o:
+
module_global_diagnostics.o:
module_OcnThickHadv.o:
@@ -182,6 +186,7 @@
                                        module_OcnVmixCoefsRich.o \
                                        module_OcnVmixCoefsTanh.o \
                                        module_OcnEquationOfState.o \
+                                        module_OcnTendency.o \
                                        module_time_integration.o
clean:
Added: branches/ocean_projects/performance/src/core_ocean/module_OcnTendency.F
===================================================================
--- branches/ocean_projects/performance/src/core_ocean/module_OcnTendency.F         (rev 0)
+++ branches/ocean_projects/performance/src/core_ocean/module_OcnTendency.F        2011-09-23 21:26:14 UTC (rev 1019)
@@ -0,0 +1,1380 @@
+!|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
+!
+! OcnTendency
+!
+!> \brief MPAS ocean tendency driver
+!> \author Doug Jacobsen
+!> \date 23 September 2011
+!> \version SVN:$Id:$
+!> \details
+!> This module contains the routines for computing
+!> various tendencies for the ocean.
+!
+!-----------------------------------------------------------------------
+
+module OcnTendency
+
+ use grid_types
+ use configure
+ use constants
+ use timer
+
+ use OcnThickHadv
+ use OcnThickVadv
+
+ use OcnVelCoriolis
+ use OcnVelPressureGrad
+ use OcnVelVadv
+ use OcnVelHmix
+ use OcnVelForcing
+
+ use OcnTracerHadv
+ use OcnTracerVadv
+ use OcnTracerHmix
+ use OcnRestoring
+
+ use OcnEquationOfState
+ use OcnVmix
+
+
+ implicit none
+ private
+ save
+
+ !--------------------------------------------------------------------
+ !
+ ! Public parameters
+ !
+ !--------------------------------------------------------------------
+
+ !--------------------------------------------------------------------
+ !
+ ! Public member functions
+ !
+ !--------------------------------------------------------------------
+
+ public :: OcnTendH, &
+ OcnTendU, &
+ OcnTendScalar, &
+ OcnDiagnosticSolve, &
+ OcnWtop, &
+ OcnFUPerp
+
+ !--------------------------------------------------------------------
+ !
+ ! Private module variables
+ !
+ !--------------------------------------------------------------------
+
+
+!***********************************************************************
+
+contains
+
+!***********************************************************************
+!
+! routine OcnTendH
+!
+!> \brief Computes thickness tendency
+!> \author Doug Jacobsen
+!> \date 23 September 2011
+!> \version SVN:$Id$
+!> \details
+!> This routine computes the thickness tendency for the ocean
+!
+!-----------------------------------------------------------------------
+
+ subroutine OcnTendH(tend, s, d, grid)!{{{
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ ! Compute height and normal wind tendencies, as well as diagnostic variables
+ !
+ ! Input: s - current model state
+ ! grid - grid metadata
+ !
+ ! Output: tend - computed tendencies for prognostic variables
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+ implicit none
+
+ type (tend_type), intent(inout) :: tend
+ type (state_type), intent(in) :: s
+ type (diagnostics_type), intent(in) :: d
+ type (mesh_type), intent(in) :: grid
+
+ integer :: iEdge, iCell, iVertex, k, cell1, cell2, &
+ vertex1, vertex2, eoe, i, j, err
+
+! mrp 110512 I just split compute_tend into compute_tend_u and OcnTendH.
+! Most of these variables can be removed, but at a later time.
+ integer :: nCells, nEdges, nVertices, nVertLevels, nEdgesSolve
+ real (kind=RKIND) :: flux, vorticity_abs, h_vertex, workpv, q, &
+ upstream_bias, wTopEdge, rho0Inv, r
+ real (kind=RKIND), dimension(:), pointer :: &
+ h_s, fVertex, fEdge, dvEdge, dcEdge, areaCell, areaTriangle, &
+ zMidZLevel, zTopZLevel
+ real (kind=RKIND), dimension(:,:), pointer :: &
+ weightsOnEdge, kiteAreasOnVertex, h_edge, h, u, v, pressure, &
+ tend_h, circulation, vorticity, ke, ke_edge, pv_edge, &
+ MontPot, wTop, divergence, vertViscTopOfEdge
+ type (dm_info) :: dminfo
+
+ integer, dimension(:), pointer :: nEdgesOnCell, nEdgesOnEdge, &
+ maxLevelCell, maxLevelEdgeTop, maxLevelVertexBot
+ integer, dimension(:,:), pointer :: &
+ cellsOnEdge, cellsOnVertex, verticesOnEdge, edgesOnCell, &
+ edgesOnEdge, edgesOnVertex
+ real (kind=RKIND) :: u_diffusion
+ real (kind=RKIND), dimension(:), allocatable:: fluxVertTop,w_dudzTopEdge
+
+ real (kind=RKIND), allocatable, dimension(:,:) :: delsq_divergence
+ real (kind=RKIND), allocatable, dimension(:,:) :: delsq_u
+ real (kind=RKIND), allocatable, dimension(:,:) :: delsq_circulation, delsq_vorticity
+
+ call timer_start("OcnTendH")
+
+ h => s % h % array
+ u => s % u % array
+ v => s % v % array
+ wTop => s % wTop % array
+ h_edge => s % h_edge % array
+ circulation => s % circulation % array
+ vorticity => s % vorticity % array
+ divergence => s % divergence % array
+ ke => s % ke % array
+ ke_edge => s % ke_edge % array
+ pv_edge => s % pv_edge % array
+ MontPot => s % MontPot % array
+ pressure => s % pressure % array
+ vertViscTopOfEdge => d % vertViscTopOfEdge % array
+
+ weightsOnEdge => grid % weightsOnEdge % array
+ kiteAreasOnVertex => grid % kiteAreasOnVertex % array
+ cellsOnEdge => grid % cellsOnEdge % array
+ cellsOnVertex => grid % cellsOnVertex % array
+ verticesOnEdge => grid % verticesOnEdge % array
+ nEdgesOnCell => grid % nEdgesOnCell % array
+ edgesOnCell => grid % edgesOnCell % array
+ nEdgesOnEdge => grid % nEdgesOnEdge % array
+ edgesOnEdge => grid % edgesOnEdge % array
+ edgesOnVertex => grid % edgesOnVertex % array
+ dcEdge => grid % dcEdge % array
+ dvEdge => grid % dvEdge % array
+ areaCell => grid % areaCell % array
+ areaTriangle => grid % areaTriangle % array
+ h_s => grid % h_s % array
+ fVertex => grid % fVertex % array
+ fEdge => grid % fEdge % array
+ zMidZLevel => grid % zMidZLevel % array
+ zTopZLevel => grid % zTopZLevel % array
+ maxLevelCell => grid % maxLevelCell % array
+ maxLevelEdgeTop => grid % maxLevelEdgeTop % array
+ maxLevelVertexBot => grid % maxLevelVertexBot % array
+
+ tend_h => tend % h % array
+
+ nCells = grid % nCells
+ nEdges = grid % nEdges
+ nEdgesSolve = grid % nEdgesSolve
+ nVertices = grid % nVertices
+ nVertLevels = grid % nVertLevels
+
+ !
+ ! height tendency: start accumulating tendency terms
+ !
+ tend_h = 0.0
+
+ !
+ ! height tendency: horizontal advection term -</font>
<font color="blue">abla\cdot ( hu)
+ !
+ ! See Ringler et al. (2010) jcp paper, eqn 19, 21, and fig. 3.
+ ! for explanation of divergence operator.
+ !
+ ! for z-level, only compute height tendency for top layer.
+
+ call timer_start("OcnTendH-horiz adv")
+
+ call OcnThickHadvTend(grid, u, h_edge, tend_h, err)
+
+ call timer_stop("OcnTendH-horiz adv")
+
+ !
+ ! height tendency: vertical advection term -d/dz(hw)
+ !
+ ! Vertical advection computed for top layer of a z grid only.
+ call timer_start("OcnTendH-vert adv")
+
+ call OcnThickVadvTend(grid, wTop, tend_h, err)
+
+ call timer_stop("OcnTendH-vert adv")
+ call timer_stop("OcnTendH")
+
+ end subroutine OcnTendH!}}}
+
+!***********************************************************************
+!
+! routine OcnTendU
+!
+!> \brief Computes velocity tendency
+!> \author Doug Jacobsen
+!> \date 23 September 2011
+!> \version SVN:$Id$
+!> \details
+!> This routine computes the velocity tendency for the ocean
+!
+!-----------------------------------------------------------------------
+
+ subroutine OcnTendU(tend, s, d, grid)!{{{
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ ! Compute height and normal wind tendencies, as well as diagnostic variables
+ !
+ ! Input: s - current model state
+ ! grid - grid metadata
+ !
+ ! Output: tend - computed tendencies for prognostic variables
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+ implicit none
+
+ type (tend_type), intent(inout) :: tend
+ type (state_type), intent(in) :: s
+ type (diagnostics_type), intent(in) :: d
+ type (mesh_type), intent(in) :: grid
+
+! mrp 110512 I just split compute_tend into OcnTendU and compute_tend_h.
+! Some of these variables can be removed, but at a later time.
+ integer :: iEdge, iCell, iVertex, k, cell1, cell2, &
+ vertex1, vertex2, eoe, i, j
+
+ integer :: nCells, nEdges, nVertices, nVertLevels, nEdgesSolve, err
+ real (kind=RKIND) :: flux, vorticity_abs, h_vertex, workpv, q, &
+ upstream_bias, wTopEdge, rho0Inv, r, visc_vorticity_coef
+ real (kind=RKIND), dimension(:), pointer :: &
+ h_s, fVertex, fEdge, dvEdge, dcEdge, areaCell, areaTriangle, &
+ zMidZLevel, zTopZLevel, meshScalingDel2, meshScalingDel4
+ real (kind=RKIND), dimension(:,:), pointer :: &
+ weightsOnEdge, kiteAreasOnVertex, h_edge, h, u, v, pressure, &
+ tend_u, circulation, vorticity, ke, ke_edge, pv_edge, &
+ MontPot, wTop, divergence, vertViscTopOfEdge
+ type (dm_info) :: dminfo
+
+ integer, dimension(:), pointer :: nEdgesOnCell, nEdgesOnEdge, &
+ maxLevelCell, maxLevelEdgeTop, maxLevelVertexBot
+ integer, dimension(:,:), pointer :: &
+ cellsOnEdge, cellsOnVertex, verticesOnEdge, edgesOnCell, &
+ edgesOnEdge, edgesOnVertex
+ real (kind=RKIND) :: u_diffusion
+ real (kind=RKIND), dimension(:), allocatable:: fluxVertTop,w_dudzTopEdge
+
+ real (kind=RKIND), allocatable, dimension(:,:) :: delsq_divergence
+ real (kind=RKIND), allocatable, dimension(:,:) :: delsq_u
+ real (kind=RKIND), allocatable, dimension(:,:) :: delsq_circulation, delsq_vorticity
+
+
+ real (kind=RKIND), dimension(:,:), pointer :: u_src
+ real (kind=RKIND), parameter :: rho_ref = 1000.0
+
+ call timer_start("OcnTendU")
+
+ h => s % h % array
+ u => s % u % array
+ v => s % v % array
+ wTop => s % wTop % array
+ h_edge => s % h_edge % array
+ circulation => s % circulation % array
+ vorticity => s % vorticity % array
+ divergence => s % divergence % array
+ ke => s % ke % array
+ ke_edge => s % ke_edge % array
+ pv_edge => s % pv_edge % array
+ MontPot => s % MontPot % array
+ pressure => s % pressure % array
+ vertViscTopOfEdge => d % vertViscTopOfEdge % array
+
+ weightsOnEdge => grid % weightsOnEdge % array
+ kiteAreasOnVertex => grid % kiteAreasOnVertex % array
+ cellsOnEdge => grid % cellsOnEdge % array
+ cellsOnVertex => grid % cellsOnVertex % array
+ verticesOnEdge => grid % verticesOnEdge % array
+ nEdgesOnCell => grid % nEdgesOnCell % array
+ edgesOnCell => grid % edgesOnCell % array
+ nEdgesOnEdge => grid % nEdgesOnEdge % array
+ edgesOnEdge => grid % edgesOnEdge % array
+ edgesOnVertex => grid % edgesOnVertex % array
+ dcEdge => grid % dcEdge % array
+ dvEdge => grid % dvEdge % array
+ areaCell => grid % areaCell % array
+ areaTriangle => grid % areaTriangle % array
+ h_s => grid % h_s % array
+! mrp 110516 cleanup fvertex fedge not used in this subroutine
+ fVertex => grid % fVertex % array
+ fEdge => grid % fEdge % array
+ zMidZLevel => grid % zMidZLevel % array
+ zTopZLevel => grid % zTopZLevel % array
+ maxLevelCell => grid % maxLevelCell % array
+ maxLevelEdgeTop => grid % maxLevelEdgeTop % array
+ maxLevelVertexBot => grid % maxLevelVertexBot % array
+
+ tend_u => tend % u % array
+
+ nCells = grid % nCells
+ nEdges = grid % nEdges
+ nEdgesSolve = grid % nEdgesSolve
+ nVertices = grid % nVertices
+ nVertLevels = grid % nVertLevels
+
+ u_src => grid % u_src % array
+
+ meshScalingDel2 => grid % meshScalingDel2 % array
+ meshScalingDel4 => grid % meshScalingDel4 % array
+
+ !
+ ! velocity tendency: start accumulating tendency terms
+ !
+ ! mrp 110516 efficiency: could remove next line and have first tend_u operation not be additive
+ tend_u(:,:) = 0.0
+
+ !
+ ! velocity tendency: nonlinear Coriolis term and grad of kinetic energy
+ !
+
+ call timer_start("OcnTendU-coriolis")
+
+! call OcnVelCoriolisTend(grid, pv_edge, h_edge, u, ke, tend_u, err)
+
+ do iEdge=1,grid % nEdgesSolve
+ cell1 = cellsOnEdge(1,iEdge)
+ cell2 = cellsOnEdge(2,iEdge)
+
+ do k=1,maxLevelEdgeTop(iEdge)
+
+ q = 0.0
+ do j = 1,nEdgesOnEdge(iEdge)
+ eoe = edgesOnEdge(j,iEdge)
+ workpv = 0.5 * (pv_edge(k,iEdge) + pv_edge(k,eoe))
+ q = q + weightsOnEdge(j,iEdge) * u(k,eoe) * workpv * h_edge(k,eoe)
+ end do
+
+ tend_u(k,iEdge) = tend_u(k,iEdge) &
+ + q &
+ - ( ke(k,cell2) - ke(k,cell1) ) / dcEdge(iEdge)
+
+ end do
+ end do
+
+ call timer_stop("OcnTendU-coriolis")
+
+ !
+ ! velocity tendency: vertical advection term -w du/dz
+ !
+ call timer_start("OcnTendU-vert adv")
+
+ call OcnVelVadvTend(grid, u, wTop, tend_u, err)
+
+ call timer_stop("OcnTendU-vert adv")
+
+ !
+ ! velocity tendency: pressure gradient
+ !
+ call timer_start("OcnTendU-pressure grad")
+
+ if (config_vert_grid_type.eq.'isopycnal') then
+ call OcnVelPressureGradTend(grid, MontPot, tend_u, err)
+ elseif (config_vert_grid_type.eq.'zlevel') then
+ call OcnVelPressureGradTend(grid, pressure, tend_u, err)
+ end if
+
+ call timer_stop("OcnTendU-pressure grad")
+
+ !
+ ! velocity tendency: del2 dissipation, </font>
<font color="black">u_2 </font>
<font color="blue">abla^2 u
+ ! computed as </font>
<font color="black">u( </font>
<font color="black">abla divergence + k \times </font>
<font color="blue">abla vorticity )
+ ! strictly only valid for config_h_mom_eddy_visc2 == constant
+ !
+ call timer_start("OcnTendU-horiz mix")
+
+ call OcnVelHmixTend(grid, divergence, vorticity, tend_u, err)
+
+ call timer_stop("OcnTendU-horiz mix")
+
+ !
+ ! velocity tendency: forcing and bottom drag
+ !
+ ! mrp 101115 note: in order to include flux boundary conditions, we will need to
+ ! know the bottom edge with nonzero velocity and place the drag there.
+
+ call timer_start("OcnTendU-forcings")
+
+ call OcnVelForcingTend(grid, u, u_src, ke_edge, h_edge, tend_u, err)
+
+ call timer_stop("OcnTendU-forcings")
+
+ !
+ ! velocity tendency: vertical mixing d/dz( nu_v du/dz))
+ !
+ if (.not.config_implicit_vertical_mix) then
+ call timer_start("OcnTendU-explicit vert mix")
+ call OcnVelVmixTendExplicit(grid, u, h_edge, vertViscTopOfEdge, tend_u, err)
+! 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
+
+! 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("OcnTendU-explicit vert mix")
+ endif
+ call timer_stop("OcnTendU")
+
+ end subroutine OcnTendU!}}}
+
+!***********************************************************************
+!
+! routine OcnTendSalar
+!
+!> \brief Computes scalar tendency
+!> \author Doug Jacobsen
+!> \date 23 September 2011
+!> \version SVN:$Id$
+!> \details
+!> This routine computes the scalar (tracer) tendency for the ocean
+!
+!-----------------------------------------------------------------------
+
+ subroutine OcnTendScalar(tend, s, d, grid)!{{{
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ !
+ ! Input: s - current model state
+ ! grid - grid metadata
+ ! note: the variable s % tracers really contains the tracers,
+ ! not tracers*h
+ !
+ ! Output: tend - computed scalar tendencies
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+ implicit none
+
+ type (tend_type), intent(inout) :: tend
+ type (state_type), intent(in) :: s
+ type (diagnostics_type), intent(in) :: d
+ type (mesh_type), intent(in) :: grid
+
+ integer :: i, k, iCell, iEdge, iTracer, cell1, cell2, upwindCell,&
+ nEdges, nCells, nCellsSolve, nVertLevels, num_tracers, err
+ real (kind=RKIND) :: invAreaCell1, invAreaCell2, tracer_turb_flux
+ real (kind=RKIND) :: flux, tracer_edge, r
+ real (kind=RKIND), dimension(:), pointer :: &
+ h_s, fVertex, fEdge, dvEdge, dcEdge, areaCell, areaTriangle
+ real (kind=RKIND), dimension(:,:), pointer :: &
+ u,h,wTop, h_edge, vertDiffTopOfCell
+ real (kind=RKIND), dimension(:,:,:), pointer :: &
+ tracers, tend_tr
+ integer, dimension(:,:), pointer :: boundaryEdge
+ type (dm_info) :: dminfo
+
+ integer, dimension(:), pointer :: nEdgesOnCell, nEdgesOnEdge, &
+ maxLevelCell, maxLevelEdgeTop, maxLevelVertexBot
+ integer, dimension(:,:), pointer :: cellsOnEdge, boundaryCell
+ real (kind=RKIND), dimension(:), pointer :: zTopZLevel,zMidZLevel, &
+ hRatioZLevelK, hRatioZLevelKm1, meshScalingDel2, meshScalingDel4
+ real (kind=RKIND), dimension(:), allocatable:: tracer2ndDer, tracersIn, tracersOut, posZMidZLevel, &
+ posZTopZLevel
+ real (kind=RKIND), dimension(:,:), allocatable:: fluxVertTop, boundaryMask
+ real (kind=RKIND), dimension(:,:,:), allocatable::tr_flux, tr_div, delsq_tracer, tracerTop
+
+
+ real (kind=RKIND) :: d2fdx2_cell1, d2fdx2_cell2
+ real (kind=RKIND), dimension(:,:,:), pointer :: deriv_two
+ real (kind=RKIND) :: coef_3rd_order, flux3Coef, cSignWTop
+
+ integer :: index_temperature, index_salinity, rrr
+ real (kind=RKIND), dimension(:), pointer :: temperatureRestore, salinityRestore
+
+ call timer_start("OcnTendScalar")
+
+ u => s % u % array
+ h => s % h % array
+ boundaryCell=> grid % boundaryCell % array
+ wTop => s % wTop % array
+ tracers => s % tracers % array
+ h_edge => s % h_edge % array
+ vertDiffTopOfCell => d % vertDiffTopOfCell % array
+
+ tend_tr => tend % tracers % array
+
+ areaCell => grid % areaCell % array
+ cellsOnEdge => grid % cellsOnEdge % array
+ dvEdge => grid % dvEdge % array
+ dcEdge => grid % dcEdge % array
+ zTopZLevel => grid % zTopZLevel % array
+ zMidZLevel => grid % zMidZLevel % array
+ hRatioZLevelK => grid % hRatioZLevelK % array
+ hRatioZLevelKm1 => grid % hRatioZLevelKm1 % array
+ boundaryEdge => grid % boundaryEdge % array
+ maxLevelCell => grid % maxLevelCell % array
+ maxLevelEdgeTop => grid % maxLevelEdgeTop % array
+ maxLevelVertexBot => grid % maxLevelVertexBot % array
+
+ nEdges = grid % nEdges
+ nCells = grid % nCells
+ nCellsSolve = grid % nCellsSolve
+ nVertLevels = grid % nVertLevels
+ num_tracers = s % num_tracers
+
+ meshScalingDel2 => grid % meshScalingDel2 % array
+ meshScalingDel4 => grid % meshScalingDel4 % array
+
+
+ deriv_two => grid % deriv_two % array
+
+ !
+ ! initialize tracer tendency (RHS of tracer equation) to zero.
+ !
+ tend_tr(:,:,:) = 0.0
+
+ !
+ ! tracer tendency: horizontal advection term -div( h \phi u)
+ !
+ ! mrp 101115 note: in order to include flux boundary conditions, we will need to
+ ! assign h_edge for maxLevelEdgeTop:maxLevelEdgeBot in the compute_solve_diagnostics
+ ! and then change maxLevelEdgeTop to maxLevelEdgeBot in the following section.
+ ! tracer_edge at the boundary will also need to be defined for flux boundaries.
+
+ call timer_start("OcnTendScalar-horiz adv")
+
+ call OcnTracerHadvTend(grid, u, h_edge, tracers, tend_tr, err)
+
+ call timer_stop("OcnTendScalar-horiz adv")
+
+
+ !
+ ! tracer tendency: vertical advection term -d/dz( h \phi w)
+ !
+
+ call timer_start("OcnTendScalar-vert adv")
+
+ call OcnTracerVadvTend(grid, wTop, tracers, tend_tr, err)
+
+ call timer_stop("OcnTendScalar-vert adv")
+
+ !
+ ! tracer tendency: del2 horizontal tracer diffusion, div(h \kappa_2 </font>
<font color="gray">abla \phi)
+ !
+ call timer_start("OcnTendScalar-horiz diff")
+
+ call OcnTracerHmixTend(grid, h_edge, tracers, tend_tr, err)
+
+ call timer_stop("OcnTendScalar-horiz diff")
+
+! mrp 110516 printing
+!print *, 'tend_tr 1',minval(tend_tr(3,1,1:nCells)),&
+! maxval(tend_tr(3,1,1:nCells))
+!print *, 'tracer 1',minval(tracers(3,1,1:nCells)),&
+! maxval(tracers(3,1,1:nCells))
+! mrp 110516 printing end
+
+ !
+ ! tracer tendency: vertical diffusion h d/dz( \kappa_v d\phi/dz)
+ !
+ if (.not.config_implicit_vertical_mix) then
+ call timer_start("OcnTendScalar-explicit vert diff")
+ call OcnTracerVmixTendExplicit(grid, h, vertDiffTopOfCell, tracers, tend_tr, err)
+! 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_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("OcnTendScalar-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))
+! mrp 110516 printing end
+
+ !
+ ! add restoring to T and S in top model layer
+ !
+ call timer_start("OcnTendScalar-restoring")
+
+ call OcnRestoringTend(grid, h, s%index_temperature, s%index_salinity, tracers, tend_tr, err)
+
+ call timer_stop("OcnTendScalar-restoring")
+
+ 10 format(2i8,10e20.10)
+ call timer_stop("OcnTendScalar")
+
+ end subroutine OcnTendScalar!}}}
+
+!***********************************************************************
+!
+! routine OcnDiagnosticSolve
+!
+!> \brief Computes diagnostic variables
+!> \author Doug Jacobsen
+!> \date 23 September 2011
+!> \version SVN:$Id$
+!> \details
+!> This routine computes the diagnostic variables for the ocean
+!
+!-----------------------------------------------------------------------
+
+ subroutine OcnDiagnosticSolve(dt, s, grid)!{{{
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ ! Compute diagnostic fields used in the tendency computations
+ !
+ ! Input: grid - grid metadata
+ !
+ ! Output: s - computed diagnostics
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+ implicit none
+
+ real (kind=RKIND), intent(in) :: dt
+ type (state_type), intent(inout) :: s
+ type (mesh_type), intent(in) :: grid
+
+
+ integer :: iEdge, iCell, iVertex, k, cell1, cell2, vertex1, vertex2, eoe, i, j, cov
+ real (kind=RKIND) :: flux, vorticity_abs, h_vertex, workpv, rho0Inv
+
+ integer :: nCells, nEdges, nVertices, nVertLevels, vertexDegree, fCoef
+
+
+ real (kind=RKIND), dimension(:), pointer :: &
+ h_s, fVertex, fEdge, dvEdge, dcEdge, areaCell, areaTriangle, &
+ hZLevel
+ real (kind=RKIND), dimension(:,:), pointer :: &
+ weightsOnEdge, kiteAreasOnVertex, h_edge, h, u, v, w, pressure,&
+ circulation, vorticity, ke, ke_edge, MontPot, wTop, &
+ pv_edge, pv_vertex, pv_cell, gradPVn, gradPVt, divergence, &
+ rho, temperature, salinity
+ real (kind=RKIND), dimension(:,:,:), pointer :: tracers
+ real (kind=RKIND), dimension(:), allocatable:: pTop
+ real (kind=RKIND), dimension(:,:), allocatable:: div_u
+ character :: c1*6
+
+ integer, dimension(:,:), pointer :: cellsOnEdge, cellsOnVertex, &
+ verticesOnEdge, edgesOnCell, edgesOnEdge, edgesOnVertex, &
+ boundaryEdge, boundaryCell
+ integer, dimension(:), pointer :: nEdgesOnCell, nEdgesOnEdge, &
+ maxLevelCell, maxLevelEdgeTop, maxLevelEdgeBot, &
+ maxLevelVertexBot, maxLevelVertexTop
+ real (kind=RKIND) :: d2fdx2_cell1, d2fdx2_cell2
+ real (kind=RKIND), dimension(:,:,:), pointer :: deriv_two
+ real (kind=RKIND) :: coef_3rd_order
+ real (kind=RKIND) :: r, h1, h2
+
+ call timer_start("OcnDiagnosticSolve")
+
+ h => s % h % array
+ u => s % u % array
+ v => s % v % array
+ wTop => s % wTop % array
+ h_edge => s % h_edge % array
+ circulation => s % circulation % array
+ vorticity => s % vorticity % array
+ divergence => s % divergence % array
+ ke => s % ke % array
+ ke_edge => s % ke_edge % array
+ pv_edge => s % pv_edge % array
+ pv_vertex => s % pv_vertex % array
+ pv_cell => s % pv_cell % array
+ gradPVn => s % gradPVn % array
+ gradPVt => s % gradPVt % array
+ rho => s % rho % array
+ tracers => s % tracers % array
+ MontPot => s % MontPot % array
+ pressure => s % pressure % array
+
+ weightsOnEdge => grid % weightsOnEdge % array
+ kiteAreasOnVertex => grid % kiteAreasOnVertex % array
+ cellsOnEdge => grid % cellsOnEdge % array
+ cellsOnVertex => grid % cellsOnVertex % array
+ verticesOnEdge => grid % verticesOnEdge % array
+ nEdgesOnCell => grid % nEdgesOnCell % array
+ edgesOnCell => grid % edgesOnCell % array
+ nEdgesOnEdge => grid % nEdgesOnEdge % array
+ edgesOnEdge => grid % edgesOnEdge % array
+ edgesOnVertex => grid % edgesOnVertex % array
+ dcEdge => grid % dcEdge % array
+ dvEdge => grid % dvEdge % array
+ areaCell => grid % areaCell % array
+ areaTriangle => grid % areaTriangle % array
+ h_s => grid % h_s % array
+ fVertex => grid % fVertex % array
+ fEdge => grid % fEdge % array
+ hZLevel => grid % hZLevel % array
+ deriv_two => grid % deriv_two % array
+ maxLevelCell => grid % maxLevelCell % array
+ maxLevelEdgeTop => grid % maxLevelEdgeTop % array
+ maxLevelEdgeBot => grid % maxLevelEdgeBot % array
+ maxLevelVertexBot => grid % maxLevelVertexBot % array
+ maxLevelVertexTop => grid % maxLevelVertexTop % array
+
+ nCells = grid % nCells
+ nEdges = grid % nEdges
+ nVertices = grid % nVertices
+ nVertLevels = grid % nVertLevels
+ vertexDegree = grid % vertexDegree
+
+ boundaryEdge => grid % boundaryEdge % array
+ boundaryCell => grid % boundaryCell % array
+
+ !
+ ! Compute height on cell edges at velocity locations
+ ! Namelist options control the order of accuracy of the reconstructed h_edge value
+ !
+ ! mrp 101115 note: in order to include flux boundary conditions, we will need to
+ ! assign h_edge for maxLevelEdgeTop:maxLevelEdgeBot in the following section
+
+ ! mrp 110516 efficiency note: For z-level, only do this on level 1. h_edge for all
+ ! lower levels is defined by hZlevel.
+
+ call timer_start("OcnDiagnosticSolve-hEdge")
+
+ coef_3rd_order = 0.
+ if (config_thickness_adv_order == 3) coef_3rd_order = 1.0
+ if (config_thickness_adv_order == 3 .and. config_monotonic) coef_3rd_order = 0.25
+
+ if (config_thickness_adv_order == 2) then
+ call timer_start("OcnDiagnosticSolve-hEdge 2")
+
+ do iEdge=1,nEdges
+ 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
+ end do
+ call timer_stop("OcnDiagnosticSolve-hEdge 2")
+
+ else if (config_thickness_adv_order == 3) then
+ call timer_start("OcnDiagnosticSolve-hEdge 3")
+
+ do iEdge=1,nEdges
+ cell1 = cellsOnEdge(1,iEdge)
+ cell2 = cellsOnEdge(2,iEdge)
+
+ do k=1,maxLevelEdgeTop(iEdge)
+
+ d2fdx2_cell1 = 0.0
+ d2fdx2_cell2 = 0.0
+
+ !-- if not a boundary cell
+ if(boundaryCell(k,cell1).eq.0.and.boundaryCell(k,cell2).eq.0) then
+
+ d2fdx2_cell1 = deriv_two(1,1,iEdge) * h(k,cell1)
+ d2fdx2_cell2 = deriv_two(1,2,iEdge) * h(k,cell2)
+
+ !-- all edges of cell 1
+ do i=1, grid % nEdgesOnCell % array (cell1)
+ d2fdx2_cell1 = d2fdx2_cell1 + &
+ deriv_two(i+1,1,iEdge) * h(k,grid % CellsOnCell % array (i,cell1))
+ end do
+
+ !-- all edges of cell 2
+ do i=1, grid % nEdgesOnCell % array (cell2)
+ d2fdx2_cell2 = d2fdx2_cell2 + &
+ deriv_two(i+1,2,iEdge) * h(k,grid % CellsOnCell % array (i,cell2))
+ end do
+
+ endif
+
+ !-- if u > 0:
+ if (u(k,iEdge) > 0) then
+ h_edge(k,iEdge) = &
+ 0.5*(h(k,cell1) + h(k,cell2)) &
+ -(dcEdge(iEdge) **2) * (d2fdx2_cell1 + d2fdx2_cell2) / 12. &
+ -(dcEdge(iEdge) **2) * coef_3rd_order*(d2fdx2_cell1 - d2fdx2_cell2) / 12.
+ !-- else u <= 0:
+ else
+ h_edge(k,iEdge) = &
+ 0.5*(h(k,cell1) + h(k,cell2)) &
+ -(dcEdge(iEdge) **2) * (d2fdx2_cell1 + d2fdx2_cell2) / 12. &
+ +(dcEdge(iEdge) **2) * coef_3rd_order*(d2fdx2_cell1 - d2fdx2_cell2) / 12.
+ end if
+
+ end do ! do k
+ end do ! do iEdge
+
+ call timer_stop("OcnDiagnosticSolve-hEdge 3")
+ else if (config_thickness_adv_order == 4) then
+ call timer_start("OcnDiagnosticSolve-hEdge 4")
+
+ do iEdge=1,nEdges
+ cell1 = cellsOnEdge(1,iEdge)
+ cell2 = cellsOnEdge(2,iEdge)
+
+ do k=1,maxLevelEdgeTop(iEdge)
+
+ d2fdx2_cell1 = 0.0
+ d2fdx2_cell2 = 0.0
+
+ !-- if not a boundary cell
+ if(boundaryCell(k,cell1).eq.0.and.boundaryCell(k,cell2).eq.0) then
+
+ d2fdx2_cell1 = deriv_two(1,1,iEdge) * h(k,cell1)
+ d2fdx2_cell2 = deriv_two(1,2,iEdge) * h(k,cell2)
+
+ !-- all edges of cell 1
+ do i=1, grid % nEdgesOnCell % array (cell1)
+ d2fdx2_cell1 = d2fdx2_cell1 + &
+ deriv_two(i+1,1,iEdge) * h(k,grid % CellsOnCell % array (i,cell1))
+ end do
+
+ !-- all edges of cell 2
+ do i=1, grid % nEdgesOnCell % array (cell2)
+ d2fdx2_cell2 = d2fdx2_cell2 + &
+ deriv_two(i+1,2,iEdge) * h(k,grid % CellsOnCell % array (i,cell2))
+ end do
+
+ endif
+
+ h_edge(k,iEdge) = &
+ 0.5*(h(k,cell1) + h(k,cell2)) &
+ -(dcEdge(iEdge) **2) * (d2fdx2_cell1 + d2fdx2_cell2) / 12.
+
+ end do ! do k
+ end do ! do iEdge
+
+ call timer_stop("OcnDiagnosticSolve-hEdge 4")
+ endif ! if(config_thickness_adv_order == 2)
+ call timer_stop("OcnDiagnosticSolve-hEdge")
+
+ !
+ ! set the velocity and height at dummy address
+ ! used -1e34 so error clearly occurs if these values are used.
+ !
+!mrp 110516 change to zero, change back later:
+ u(:,nEdges+1) = -1e34
+ h(:,nCells+1) = -1e34
+ tracers(s % index_temperature,:,nCells+1) = -1e34
+ tracers(s % index_salinity,:,nCells+1) = -1e34
+
+ !
+ ! Compute circulation and relative vorticity at each vertex
+ !
+ circulation(:,:) = 0.0
+ do iEdge=1,nEdges
+ vertex1 = verticesOnEdge(1,iEdge)
+ vertex2 = verticesOnEdge(2,iEdge)
+ do k=1,maxLevelEdgeBot(iEdge)
+ circulation(k,vertex1) = circulation(k,vertex1) - dcEdge(iEdge) * u(k,iEdge)
+ circulation(k,vertex2) = circulation(k,vertex2) + dcEdge(iEdge) * u(k,iEdge)
+ end do
+ end do
+ do iVertex=1,nVertices
+ do k=1,maxLevelVertexBot(iVertex)
+ vorticity(k,iVertex) = circulation(k,iVertex) / areaTriangle(iVertex)
+ end do
+ end do
+
+ !
+ ! Compute the divergence at each cell center
+ !
+ divergence(:,:) = 0.0
+ do iEdge=1,nEdges
+ cell1 = cellsOnEdge(1,iEdge)
+ cell2 = cellsOnEdge(2,iEdge)
+ do k=1,maxLevelEdgeBot(iEdge)
+ divergence(k,cell1) = divergence(k,cell1) + u(k,iEdge)*dvEdge(iEdge)
+ divergence(k,cell2) = divergence(k,cell2) - u(k,iEdge)*dvEdge(iEdge)
+ enddo
+ end do
+ do iCell = 1,nCells
+ r = 1.0 / areaCell(iCell)
+ do k = 1,maxLevelCell(iCell)
+ divergence(k,iCell) = divergence(k,iCell) * r
+ enddo
+ enddo
+
+ !
+ ! Compute kinetic energy in each cell
+ !
+ ke(:,:) = 0.0
+ do iEdge=1,nEdges
+ cell1 = cellsOnEdge(1,iEdge)
+ cell2 = cellsOnEdge(2,iEdge)
+ do k=1,maxLevelEdgeBot(iEdge)
+ ke(k,cell1) = ke(k,cell1) + 0.25 * dcEdge(iEdge) * dvEdge(iEdge) * u(k,iEdge)**2.0
+ ke(k,cell2) = ke(k,cell2) + 0.25 * dcEdge(iEdge) * dvEdge(iEdge) * u(k,iEdge)**2.0
+ enddo
+ end do
+ do iCell = 1,nCells
+ do k = 1,maxLevelCell(iCell)
+ ke(k,iCell) = ke(k,iCell) / areaCell(iCell)
+ enddo
+ enddo
+
+ !
+ ! Compute v (tangential) velocities
+ !
+ v(:,:) = 0.0
+ do iEdge = 1,nEdges
+ do i=1,nEdgesOnEdge(iEdge)
+ eoe = edgesOnEdge(i,iEdge)
+ ! mrp 101115 note: in order to include flux boundary conditions,
+ ! the following loop may need to change to maxLevelEdgeBot
+ do k = 1,maxLevelEdgeTop(iEdge)
+ v(k,iEdge) = v(k,iEdge) + weightsOnEdge(i,iEdge) * u(k, eoe)
+ end do
+ end do
+ end do
+
+ !
+ ! Compute ke on cell edges at velocity locations for quadratic bottom drag.
+ !
+ ! mrp 101025 efficiency note: we could get rid of ke_edge completely by
+ ! using sqrt(u(k,iEdge)**2 + v(k,iEdge)**2) in its place elsewhere.
+ ke_edge = 0.0 !mrp remove 0 for efficiency
+ do iEdge=1,nEdges
+ cell1 = cellsOnEdge(1,iEdge)
+ cell2 = cellsOnEdge(2,iEdge)
+ do k=1,maxLevelEdgeTop(iEdge)
+ ke_edge(k,iEdge) = 0.5 * (ke(k,cell1) + ke(k,cell2))
+ end do
+ end do
+
+ !
+ ! Compute height at vertices, pv at vertices, and average pv to edge locations
+ ! ( this computes pv_vertex at all vertices bounding real cells and distance-1 ghost cells )
+ !
+ if (trim(config_time_integration) == 'RK4') then
+ ! for RK4, PV is really PV = (eta+f)/h
+ fCoef = 1
+ elseif (trim(config_time_integration) == 'split_explicit' &
+ .or.trim(config_time_integration) == 'unsplit_explicit') then
+ ! for split explicit, PV is eta/h because f is added separately to the momentum forcing.
+! mrp temp, new should be:
+ fCoef = 0
+! old, for testing:
+! fCoef = 1
+ end if
+
+ do iVertex = 1,nVertices
+ do k=1,maxLevelVertexBot(iVertex)
+ h_vertex = 0.0
+ do i=1,vertexDegree
+ h_vertex = h_vertex + h(k,cellsOnVertex(i,iVertex)) * kiteAreasOnVertex(i,iVertex)
+ end do
+ h_vertex = h_vertex / areaTriangle(iVertex)
+
+ pv_vertex(k,iVertex) = (fCoef*fVertex(iVertex) + vorticity(k,iVertex)) / h_vertex
+ end do
+ end do
+
+ !
+ ! Compute pv at cell centers
+ ! ( this computes pv_cell for all real cells and distance-1 ghost cells )
+ !
+ pv_cell(:,:) = 0.0
+ do iVertex = 1,nVertices
+ do i=1,vertexDegree
+ iCell = cellsOnVertex(i,iVertex)
+ do k = 1,maxLevelCell(iCell)
+ pv_cell(k,iCell) = pv_cell(k,iCell) &
+ + kiteAreasOnVertex(i, iVertex) * pv_vertex(k, iVertex) &
+ / areaCell(iCell)
+ enddo
+ enddo
+ enddo
+
+ !
+ ! Compute pv at the edges
+ ! ( this computes pv_edge at all edges bounding real cells )
+ !
+ pv_edge(:,:) = 0.0
+ do iVertex = 1,nVertices
+ do i=1,vertexDegree
+ iEdge = edgesOnVertex(i,iVertex)
+ do k=1,maxLevelEdgeBot(iEdge)
+ pv_edge(k,iEdge) = pv_edge(k,iEdge) + 0.5 * pv_vertex(k,iVertex)
+ enddo
+ end do
+ end do
+
+ !
+ ! Compute gradient of PV in normal direction
+ ! ( this computes gradPVn for all edges bounding real cells )
+ !
+ gradPVn(:,:) = 0.0
+ do iEdge = 1,nEdges
+ do k=1,maxLevelEdgeTop(iEdge)
+ gradPVn(k,iEdge) = ( pv_cell(k,cellsOnEdge(2,iEdge)) &
+ - pv_cell(k,cellsOnEdge(1,iEdge))) &
+ / dcEdge(iEdge)
+ enddo
+ enddo
+
+ !
+ ! Compute gradient of PV in the tangent direction
+ ! ( this computes gradPVt at all edges bounding real cells and distance-1 ghost cells )
+ !
+ do iEdge = 1,nEdges
+ do k = 1,maxLevelEdgeBot(iEdge)
+ gradPVt(k,iEdge) = ( pv_vertex(k,verticesOnEdge(2,iEdge)) &
+ - pv_vertex(k,verticesOnEdge(1,iEdge))) &
+ /dvEdge(iEdge)
+ enddo
+ enddo
+
+ !
+ ! Modify PV edge with upstream bias.
+ !
+ do iEdge = 1,nEdges
+ do k = 1,maxLevelEdgeBot(iEdge)
+ pv_edge(k,iEdge) = pv_edge(k,iEdge) &
+ - 0.5 * dt* ( u(k,iEdge) * gradPVn(k,iEdge) &
+ + v(k,iEdge) * gradPVt(k,iEdge) )
+ enddo
+ enddo
+
+ !
+ ! equation of state
+ !
+ ! For an isopycnal model, density should remain constant.
+ ! For zlevel, calculate in-situ density
+ 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
+ !
+ if (config_vert_grid_type.eq.'isopycnal') then
+
+ ! For Isopycnal model.
+ ! Compute pressure at top of each layer, and then
+ ! Montgomery Potential.
+ allocate(pTop(nVertLevels))
+ do iCell=1,nCells
+
+ ! assume atmospheric pressure at the surface is zero for now.
+ pTop(1) = 0.0
+ ! For isopycnal mode, p is the Montgomery Potential.
+ ! At top layer it is g*SSH, where SSH may be off by a
+ ! constant (ie, h_s can be relative to top or bottom)
+ MontPot(1,iCell) = gravity &
+ * (h_s(iCell) + sum(h(1:nVertLevels,iCell)))
+
+ do k=2,nVertLevels
+ pTop(k) = pTop(k-1) + rho(k-1,iCell)*gravity* h(k-1,iCell)
+
+ ! from delta M = p delta / rho
+ MontPot(k,iCell) = MontPot(k-1,iCell) &
+ + pTop(k)*(1.0/rho(k,iCell) - 1.0/rho(k-1,iCell))
+ end do
+
+ end do
+ deallocate(pTop)
+
+ elseif (config_vert_grid_type.eq.'zlevel') then
+
+ ! For z-level model.
+ ! Compute pressure at middle of each level.
+ ! At k=1, where p is pressure at a depth of hZLevel(1)/2, not
+ ! pressure at middle of layer including SSH.
+
+ do iCell=1,nCells
+ ! compute pressure for z-level coordinates
+ ! assume atmospheric pressure at the surface is zero for now.
+
+ pressure(1,iCell) = rho(1,iCell)*gravity &
+ * (h(1,iCell)-0.5*hZLevel(1))
+
+ do k=2,maxLevelCell(iCell)
+ pressure(k,iCell) = pressure(k-1,iCell) &
+ + 0.5*gravity*( rho(k-1,iCell)*hZLevel(k-1) &
+ + rho(k ,iCell)*hZLevel(k ))
+ end do
+
+ end do
+
+ endif
+
+ call OcnWtop(s,grid)
+
+ call timer_stop("OcnDiagnosticSolve")
+
+ end subroutine OcnDiagnosticSolve!}}}
+
+!***********************************************************************
+!
+! routine OcnWtop
+!
+!> \brief Computes vertical velocity
+!> \author Doug Jacobsen
+!> \date 23 September 2011
+!> \version SVN:$Id$
+!> \details
+!> This routine computes the vertical velocity in the top layer for the ocean
+!
+!-----------------------------------------------------------------------
+
+ subroutine OcnWtop(s, grid)!{{{
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ ! Compute diagnostic fields used in the tendency computations
+ !
+ ! Input: grid - grid metadata
+ !
+ ! Output: s - computed diagnostics
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+ implicit none
+
+ type (state_type), intent(inout) :: s
+ type (mesh_type), intent(in) :: grid
+
+ ! mrp 110512 could clean this out, remove pointers?
+ integer :: iEdge, iCell, iVertex, k, cell1, cell2, vertex1, vertex2, eoe, i, j, cov
+ real (kind=RKIND) :: flux, vorticity_abs, h_vertex, workpv, rho0Inv
+
+ integer :: nCells, nEdges, nVertices, nVertLevels, vertexDegree
+
+
+ real (kind=RKIND), dimension(:), pointer :: &
+ h_s, fVertex, fEdge, dvEdge, dcEdge, areaCell, areaTriangle, &
+ hZLevel
+ real (kind=RKIND), dimension(:,:), pointer :: u,wTop
+ real (kind=RKIND), dimension(:,:), allocatable:: div_u
+
+ integer, dimension(:,:), pointer :: cellsOnEdge, cellsOnVertex, &
+ verticesOnEdge, edgesOnCell, edgesOnEdge, edgesOnVertex, &
+ boundaryEdge, boundaryCell
+ integer, dimension(:), pointer :: nEdgesOnCell, nEdgesOnEdge, &
+ maxLevelCell, maxLevelEdgeTop, maxLevelEdgeBot, &
+ maxLevelVertexBot, maxLevelVertexTop
+
+ call timer_start("wTop")
+
+ u => s % u % array
+ wTop => s % wTop % array
+
+ areaCell => grid % areaCell % array
+ cellsOnEdge => grid % cellsOnEdge % array
+ hZLevel => grid % hZLevel % array
+ maxLevelCell => grid % maxLevelCell % array
+ maxLevelEdgeBot => grid % maxLevelEdgeBot % array
+ dvEdge => grid % dvEdge % array
+
+ nCells = grid % nCells
+ nEdges = grid % nEdges
+ nVertLevels = grid % nVertLevels
+
+ !
+ ! vertical velocity through layer interface
+ !
+ if (config_vert_grid_type.eq.'isopycnal') then
+ ! set vertical velocity to zero in isopycnal case
+ wTop=0.0
+
+ elseif (config_vert_grid_type.eq.'zlevel') then
+
+ !
+ ! Compute div(u) for each cell
+ ! See Ringler et al. (2010) jcp paper, eqn 19, 21, and fig. 3.
+ !
+ allocate(div_u(nVertLevels,nCells+1))
+ div_u(:,:) = 0.0
+ do iEdge=1,nEdges
+ cell1 = cellsOnEdge(1,iEdge)
+ cell2 = cellsOnEdge(2,iEdge)
+ do k=2,maxLevelEdgeBot(iEdge)
+ flux = u(k,iEdge) * dvEdge(iEdge)
+ div_u(k,cell1) = div_u(k,cell1) + flux
+ div_u(k,cell2) = div_u(k,cell2) - flux
+ end do
+ end do
+
+ do iCell=1,nCells
+ ! Vertical velocity through layer interface at top and
+ ! bottom is zero.
+ wTop(1,iCell) = 0.0
+ wTop(maxLevelCell(iCell)+1,iCell) = 0.0
+ do k=maxLevelCell(iCell),2,-1
+ wTop(k,iCell) = wTop(k+1,iCell) &
+ - div_u(k,iCell)/areaCell(iCell)*hZLevel(k)
+ end do
+ end do
+ deallocate(div_u)
+
+ endif
+
+ call timer_stop("wTop")
+
+ end subroutine OcnWtop!}}}
+
+!***********************************************************************
+!
+! routine OcnFUPerp
+!
+!> \brief Computes f u_perp
+!> \author Doug Jacobsen
+!> \date 23 September 2011
+!> \version SVN:$Id$
+!> \details
+!> This routine computes f u_perp for the ocean
+!
+!-----------------------------------------------------------------------
+
+ subroutine OcnFUPerp(s, grid)!{{{
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ ! Put f*uBcl^{perp} in u as a work variable
+ !
+ ! Input: s - current model state
+ ! grid - grid metadata
+ !
+ ! Output: tend - computed tendencies for prognostic variables
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+ implicit none
+
+ type (state_type), intent(inout) :: s
+ type (mesh_type), intent(in) :: grid
+
+! mrp 110512 I just split compute_tend into compute_tend_u and compute_tend_h.
+! Some of these variables can be removed, but at a later time.
+ integer :: iEdge, iCell, iVertex, k, cell1, cell2, &
+ vertex1, vertex2, eoe, i, j
+
+ integer :: nCells, nEdges, nVertices, nVertLevels, nEdgesSolve
+ real (kind=RKIND) :: flux, vorticity_abs, h_vertex, workpv, q, &
+ upstream_bias, wTopEdge, rho0Inv, r
+ real (kind=RKIND), dimension(:), pointer :: &
+ h_s, fVertex, fEdge, dvEdge, dcEdge, areaCell, areaTriangle, &
+ zMidZLevel, zTopZLevel
+ real (kind=RKIND), dimension(:,:), pointer :: &
+ weightsOnEdge, kiteAreasOnVertex, h_edge, h, u, uBcl, v, pressure, &
+ tend_u, circulation, vorticity, ke, ke_edge, pv_edge, &
+ MontPot, wTop, divergence, vertViscTopOfEdge
+ type (dm_info) :: dminfo
+
+ integer, dimension(:), pointer :: nEdgesOnCell, nEdgesOnEdge, &
+ maxLevelCell, maxLevelEdgeTop, maxLevelVertexBot
+ integer, dimension(:,:), pointer :: &
+ cellsOnEdge, cellsOnVertex, verticesOnEdge, edgesOnCell, &
+ edgesOnEdge, edgesOnVertex
+ real (kind=RKIND) :: u_diffusion
+ real (kind=RKIND), dimension(:), allocatable:: fluxVertTop,w_dudzTopEdge
+
+ real (kind=RKIND), allocatable, dimension(:,:) :: delsq_divergence
+ real (kind=RKIND), allocatable, dimension(:,:) :: delsq_u
+ real (kind=RKIND), allocatable, dimension(:,:) :: delsq_circulation, delsq_vorticity
+
+
+ real (kind=RKIND), dimension(:,:), pointer :: u_src
+ real (kind=RKIND), parameter :: rho_ref = 1000.0
+
+ call timer_start("OcnFUPerp")
+
+ h => s % h % array
+ u => s % u % array
+ uBcl => s % uBcl % array
+ v => s % v % array
+ wTop => s % wTop % array
+ h_edge => s % h_edge % array
+ circulation => s % circulation % array
+ vorticity => s % vorticity % array
+ divergence => s % divergence % array
+ ke => s % ke % array
+ ke_edge => s % ke_edge % array
+ pv_edge => s % pv_edge % array
+ MontPot => s % MontPot % array
+ pressure => s % pressure % array
+
+ weightsOnEdge => grid % weightsOnEdge % array
+ kiteAreasOnVertex => grid % kiteAreasOnVertex % array
+ cellsOnEdge => grid % cellsOnEdge % array
+ cellsOnVertex => grid % cellsOnVertex % array
+ verticesOnEdge => grid % verticesOnEdge % array
+ nEdgesOnCell => grid % nEdgesOnCell % array
+ edgesOnCell => grid % edgesOnCell % array
+ nEdgesOnEdge => grid % nEdgesOnEdge % array
+ edgesOnEdge => grid % edgesOnEdge % array
+ edgesOnVertex => grid % edgesOnVertex % array
+ dcEdge => grid % dcEdge % array
+ dvEdge => grid % dvEdge % array
+ areaCell => grid % areaCell % array
+ areaTriangle => grid % areaTriangle % array
+ h_s => grid % h_s % array
+ fVertex => grid % fVertex % array
+ fEdge => grid % fEdge % array
+ zMidZLevel => grid % zMidZLevel % array
+ zTopZLevel => grid % zTopZLevel % array
+ maxLevelCell => grid % maxLevelCell % array
+ maxLevelEdgeTop => grid % maxLevelEdgeTop % array
+ maxLevelVertexBot => grid % maxLevelVertexBot % array
+
+
+ nCells = grid % nCells
+ nEdges = grid % nEdges
+ nEdgesSolve = grid % nEdgesSolve
+ nVertices = grid % nVertices
+ nVertLevels = grid % nVertLevels
+
+ !
+ ! Put f*uBcl^{perp} in u as a work variable
+ !
+ do iEdge=1,grid % nEdgesSolve
+ cell1 = cellsOnEdge(1,iEdge)
+ cell2 = cellsOnEdge(2,iEdge)
+
+ do k=1,maxLevelEdgeTop(iEdge)
+
+ u(k,iEdge) = 0.0
+ do j = 1,nEdgesOnEdge(iEdge)
+ eoe = edgesOnEdge(j,iEdge)
+ u(k,iEdge) = u(k,iEdge) + weightsOnEdge(j,iEdge) * uBcl(k,eoe) * fEdge(eoe)
+ end do
+ end do
+ end do
+
+ call timer_stop("OcnFUPerp")
+
+ end subroutine OcnFUPerp!}}}
+
+!***********************************************************************
+
+end module OcnTendency
+
+!|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
+! 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-23 18:09:34 UTC (rev 1018)
+++ branches/ocean_projects/performance/src/core_ocean/module_mpas_core.F        2011-09-23 21:26:14 UTC (rev 1019)
@@ -5,6 +5,8 @@
use dmpar
use test_cases
+ use OcnTendency
+
use OcnVelPressureGrad
use OcnVelVadv
use OcnVelHmix
@@ -188,7 +190,7 @@
integer :: i, iEdge, iCell, k
- call compute_solve_diagnostics(dt, block % state % time_levs(1) % state, mesh)
+ call OcnDiagnosticSolve(dt, block % state % time_levs(1) % state, mesh)
call compute_mesh_scaling(mesh)
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-23 18:09:34 UTC (rev 1018)
+++ branches/ocean_projects/performance/src/core_ocean/module_time_integration.F        2011-09-23 21:26:14 UTC (rev 1019)
@@ -8,20 +8,8 @@
use spline_interpolation
use timer
- use OcnThickHadv
- use OcnThickVadv
+ use OcnTendency
- use OcnVelCoriolis
- use OcnVelPressureGrad
- use OcnVelVadv
- use OcnVelHmix
- use OcnVelForcing
-
- use OcnTracerHadv
- use OcnTracerVadv
- use OcnTracerHmix
- use OcnRestoring
-
use OcnEquationOfState
use OcnVmix
@@ -182,8 +170,8 @@
if (.not.config_implicit_vertical_mix) then
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)
+ call OcnTendH(block % tend, provis, block % diagnostics, block % mesh)
+ call OcnTendU(block % tend, 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
@@ -191,7 +179,7 @@
call filter_btr_mode_tend_u(block % tend, provis, block % diagnostics, block % mesh)
endif
- call compute_scalar_tend(block % tend, provis, block % diagnostics, block % mesh)
+ call OcnTendScalar(block % tend, provis, block % diagnostics, block % mesh)
call enforce_boundaryEdge(block % tend, block % mesh)
block => block % next
end do
@@ -241,7 +229,7 @@
provis % u % array(:,:) = block % state % time_levs(1) % state % u % array(:,:)
end if
- call compute_solve_diagnostics(dt, provis, block % mesh)
+ call OcnDiagnosticSolve(dt, provis, block % mesh)
block => block % next
end do
@@ -318,41 +306,7 @@
! Implicit vertical solve for momentum
!
call OcnVelVmixTendImplicit(block % mesh, dt, ke_edge, vertViscTopOfEdge, h, h_edge, u, err)
-! 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
-
! 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)
@@ -364,32 +318,6 @@
!
call OcnTracerVmixTendImplicit(block % mesh, dt, vertDiffTopOfCell, h, tracers, err)
-! 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")
end if
! mrp 110725 momentum decay term
@@ -419,7 +347,7 @@
block % state % time_levs(2) % state % u % array(:,:) = block % state % time_levs(1) % state % u % array(:,:)
end if
- call compute_solve_diagnostics(dt, block % state % time_levs(2) % state, block % mesh)
+ call OcnDiagnosticSolve(dt, block % state % time_levs(2) % state, block % mesh)
call reconstruct(block % state % time_levs(2) % state, block % mesh)
@@ -557,7 +485,7 @@
if (.not.config_implicit_vertical_mix) then
call OcnVmixCoefs(block % mesh, block % state % time_levs(2) % state, block % diagnostics, err)
end if
- call compute_tend_u(block % tend, block % state % time_levs(2) % state , block % diagnostics, block % mesh)
+ call OcnTendU(block % tend, block % state % time_levs(2) % state , block % diagnostics, block % mesh)
call enforce_boundaryEdge(block % tend, block % mesh)
block => block % next
end do
@@ -579,7 +507,7 @@
allocate(uTemp(block % mesh % nVertLevels))
! Put f*uBcl^{perp} in uNew as a work variable
- call compute_f_uperp(block % state % time_levs(2) % state , block % mesh)
+ call OcnFUPerp(block % state % time_levs(2) % state , block % mesh)
do iEdge=1,block % mesh % nEdges
cell1 = block % mesh % cellsOnEdge % array(1,iEdge)
@@ -1189,7 +1117,7 @@
block % state % time_levs(2) % state % u % array(:,iEdge) = uTemp(:) + uCorr
endif
- ! Put new sshEdge values in h_edge array, for the compute_scalar_tend call below.
+ ! Put new sshEdge values in h_edge array, for the OcnTendScalar call below.
block % state % time_levs(2) % state % h_edge % array(1,iEdge) &
= sshEdge + block % mesh % hZLevel % array(1)
@@ -1200,7 +1128,7 @@
end do ! iEdge
- ! Put new SSH values in h array, for the compute_scalar_tend call below.
+ ! Put new SSH values in h array, for the OcnTendScalar call below.
do iCell=1,block % mesh % nCells
block % state % time_levs(2) % state % h % array(1,iCell) &
= block % state % time_levs(2) % state % ssh % array(iCell) &
@@ -1231,13 +1159,13 @@
block => domain % blocklist
do while (associated(block))
- call compute_wTop(block % state % time_levs(2) % state, block % mesh)
+ call OcnWtop(block % state % time_levs(2) % state, block % mesh)
if (trim(config_time_integration) == 'unsplit_explicit') then
- call compute_tend_h(block % tend, block % state % time_levs(2) % state , block % diagnostics, block % mesh)
+ call OcnTendH(block % tend, block % state % time_levs(2) % state , block % diagnostics, block % mesh)
endif
- call compute_scalar_tend(block % tend, block % state % time_levs(2) % state , block % diagnostics, block % mesh)
+ call OcnTendScalar(block % tend, block % state % time_levs(2) % state , block % diagnostics, block % mesh)
block => block % next
end do
@@ -1328,7 +1256,7 @@
! compute u*, the velocity for tendency terms. Put in uNew.
! uBclNew is at time n+1/2 here.
! This overwrites u^{tr}, the tracer transport velocity, which was in uNew.
- ! The following must occur after call compute_scalar_tend
+ ! The following must occur after call OcnTendScalar
do iEdge=1,block % mesh % nEdges
block % state % time_levs(2) % state % u % array(:,iEdge) &
= block % state % time_levs(2) % state % uBtr % array(iEdge) &
@@ -1337,7 +1265,7 @@
! mrp 110512 I really only need this to compute h_edge, density, pressure.
! I can par this down later.
- call compute_solve_diagnostics(dt, block % state % time_levs(2) % state, block % mesh)
+ call OcnDiagnosticSolve(dt, block % state % time_levs(2) % state, block % mesh)
elseif (split_explicit_step == config_n_ts_iter) then
@@ -1430,7 +1358,7 @@
endif
do iCell=1,block % mesh % nCells
- ! Put new SSH values in h array, for the compute_scalar_tend call below.
+ ! Put new SSH values in h array, for the OcnTendScalar call below.
block % state % time_levs(2) % state % h % array(1,iCell) &
= block % state % time_levs(2) % state % ssh % array(iCell) &
+ block % mesh % hZLevel % array(1)
@@ -1473,69 +1401,10 @@
call OcnVelVmixTendImplicit(block % mesh, dt, ke_edge, vertViscTopOfEdge, h, h_edge, u, err)
-! do iEdge=1,block % mesh % 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)) = 0.0
-
-! 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:block % mesh % nVertLevels,iEdge) = 0.0
-
-! end if
-! end do
-
!
! Implicit vertical solve for tracers
!
call OcnTracerVmixTendImplicit(block % mesh, dt, vertDiffTopOfCell, h, tracers, err)
-! do iCell=1,block % mesh % 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), &
-! block % mesh % nVertLevels,num_tracers)
-
-! tracers(:,1:maxLevelCell(iCell),iCell) = tracersTemp(:,1:maxLevelCell(iCell))
-! tracers(:,maxLevelCell(iCell)+1:block % mesh % nVertLevels,iCell) = -1e34
-! end do
-! deallocate(A,C,uTemp,tracersTemp)
end if
! mrp 110725 adding momentum decay term
@@ -1562,7 +1431,7 @@
block % state % time_levs(2) % state % u % array(:,:) = block % state % time_levs(1) % state % u % array(:,:)
end if
- call compute_solve_diagnostics(dt, block % state % time_levs(2) % state, block % mesh)
+ call OcnDiagnosticSolve(dt, block % state % time_levs(2) % state, block % mesh)
call reconstruct(block % state % time_levs(2) % state, block % mesh)
@@ -1572,348 +1441,6 @@
end subroutine split_explicit_timestep!}}}
- subroutine compute_tend_h(tend, s, d, grid)!{{{
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
- ! Compute height and normal wind tendencies, as well as diagnostic variables
- !
- ! Input: s - current model state
- ! grid - grid metadata
- !
- ! Output: tend - computed tendencies for prognostic variables
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
- implicit none
-
- type (tend_type), intent(inout) :: tend
- type (state_type), intent(in) :: s
- type (diagnostics_type), intent(in) :: d
- type (mesh_type), intent(in) :: grid
-
- integer :: iEdge, iCell, iVertex, k, cell1, cell2, &
- vertex1, vertex2, eoe, i, j, err
-
-! mrp 110512 I just split compute_tend into compute_tend_u and compute_tend_h.
-! Most of these variables can be removed, but at a later time.
- integer :: nCells, nEdges, nVertices, nVertLevels, nEdgesSolve
- real (kind=RKIND) :: flux, vorticity_abs, h_vertex, workpv, q, &
- upstream_bias, wTopEdge, rho0Inv, r
- real (kind=RKIND), dimension(:), pointer :: &
- h_s, fVertex, fEdge, dvEdge, dcEdge, areaCell, areaTriangle, &
- zMidZLevel, zTopZLevel
- real (kind=RKIND), dimension(:,:), pointer :: &
- weightsOnEdge, kiteAreasOnVertex, h_edge, h, u, v, pressure, &
- tend_h, circulation, vorticity, ke, ke_edge, pv_edge, &
- MontPot, wTop, divergence, vertViscTopOfEdge
- type (dm_info) :: dminfo
-
- integer, dimension(:), pointer :: nEdgesOnCell, nEdgesOnEdge, &
- maxLevelCell, maxLevelEdgeTop, maxLevelVertexBot
- integer, dimension(:,:), pointer :: &
- cellsOnEdge, cellsOnVertex, verticesOnEdge, edgesOnCell, &
- edgesOnEdge, edgesOnVertex
- real (kind=RKIND) :: u_diffusion
- real (kind=RKIND), dimension(:), allocatable:: fluxVertTop,w_dudzTopEdge
-
- real (kind=RKIND), allocatable, dimension(:,:) :: delsq_divergence
- real (kind=RKIND), allocatable, dimension(:,:) :: delsq_u
- real (kind=RKIND), allocatable, dimension(:,:) :: delsq_circulation, delsq_vorticity
-
- call timer_start("compute_tend_h")
-
- h => s % h % array
- u => s % u % array
- v => s % v % array
- wTop => s % wTop % array
- h_edge => s % h_edge % array
- circulation => s % circulation % array
- vorticity => s % vorticity % array
- divergence => s % divergence % array
- ke => s % ke % array
- ke_edge => s % ke_edge % array
- pv_edge => s % pv_edge % array
- MontPot => s % MontPot % array
- pressure => s % pressure % array
- vertViscTopOfEdge => d % vertViscTopOfEdge % array
-
- weightsOnEdge => grid % weightsOnEdge % array
- kiteAreasOnVertex => grid % kiteAreasOnVertex % array
- cellsOnEdge => grid % cellsOnEdge % array
- cellsOnVertex => grid % cellsOnVertex % array
- verticesOnEdge => grid % verticesOnEdge % array
- nEdgesOnCell => grid % nEdgesOnCell % array
- edgesOnCell => grid % edgesOnCell % array
- nEdgesOnEdge => grid % nEdgesOnEdge % array
- edgesOnEdge => grid % edgesOnEdge % array
- edgesOnVertex => grid % edgesOnVertex % array
- dcEdge => grid % dcEdge % array
- dvEdge => grid % dvEdge % array
- areaCell => grid % areaCell % array
- areaTriangle => grid % areaTriangle % array
- h_s => grid % h_s % array
- fVertex => grid % fVertex % array
- fEdge => grid % fEdge % array
- zMidZLevel => grid % zMidZLevel % array
- zTopZLevel => grid % zTopZLevel % array
- maxLevelCell => grid % maxLevelCell % array
- maxLevelEdgeTop => grid % maxLevelEdgeTop % array
- maxLevelVertexBot => grid % maxLevelVertexBot % array
-
- tend_h => tend % h % array
-
- nCells = grid % nCells
- nEdges = grid % nEdges
- nEdgesSolve = grid % nEdgesSolve
- nVertices = grid % nVertices
- nVertLevels = grid % nVertLevels
-
- !
- ! height tendency: start accumulating tendency terms
- !
- tend_h = 0.0
-
- !
- ! height tendency: horizontal advection term -</font>
<font color="red">abla\cdot ( hu)
- !
- ! See Ringler et al. (2010) jcp paper, eqn 19, 21, and fig. 3.
- ! for explanation of divergence operator.
- !
- ! for z-level, only compute height tendency for top layer.
-
- call timer_start("compute_tend_h-horiz adv")
-
- call OcnThickHadvTend(grid, u, h_edge, tend_h, err)
-
- call timer_stop("compute_tend_h-horiz adv")
-
- !
- ! height tendency: vertical advection term -d/dz(hw)
- !
- ! Vertical advection computed for top layer of a z grid only.
- call timer_start("compute_tend_h-vert adv")
-
- call OcnThickVadvTend(grid, wTop, tend_h, err)
-
- call timer_stop("compute_tend_h-vert adv")
- call timer_stop("compute_tend_h")
-
- end subroutine compute_tend_h!}}}
-
- subroutine compute_tend_u(tend, s, d, grid)!{{{
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
- ! Compute height and normal wind tendencies, as well as diagnostic variables
- !
- ! Input: s - current model state
- ! grid - grid metadata
- !
- ! Output: tend - computed tendencies for prognostic variables
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
- implicit none
-
- type (tend_type), intent(inout) :: tend
- type (state_type), intent(in) :: s
- type (diagnostics_type), intent(in) :: d
- type (mesh_type), intent(in) :: grid
-
-! mrp 110512 I just split compute_tend into compute_tend_u and compute_tend_h.
-! Some of these variables can be removed, but at a later time.
- integer :: iEdge, iCell, iVertex, k, cell1, cell2, &
- vertex1, vertex2, eoe, i, j
-
- integer :: nCells, nEdges, nVertices, nVertLevels, nEdgesSolve, err
- real (kind=RKIND) :: flux, vorticity_abs, h_vertex, workpv, q, &
- upstream_bias, wTopEdge, rho0Inv, r, visc_vorticity_coef
- real (kind=RKIND), dimension(:), pointer :: &
- h_s, fVertex, fEdge, dvEdge, dcEdge, areaCell, areaTriangle, &
- zMidZLevel, zTopZLevel, meshScalingDel2, meshScalingDel4
- real (kind=RKIND), dimension(:,:), pointer :: &
- weightsOnEdge, kiteAreasOnVertex, h_edge, h, u, v, pressure, &
- tend_u, circulation, vorticity, ke, ke_edge, pv_edge, &
- MontPot, wTop, divergence, vertViscTopOfEdge
- type (dm_info) :: dminfo
-
- integer, dimension(:), pointer :: nEdgesOnCell, nEdgesOnEdge, &
- maxLevelCell, maxLevelEdgeTop, maxLevelVertexBot
- integer, dimension(:,:), pointer :: &
- cellsOnEdge, cellsOnVertex, verticesOnEdge, edgesOnCell, &
- edgesOnEdge, edgesOnVertex
- real (kind=RKIND) :: u_diffusion
- real (kind=RKIND), dimension(:), allocatable:: fluxVertTop,w_dudzTopEdge
-
- real (kind=RKIND), allocatable, dimension(:,:) :: delsq_divergence
- real (kind=RKIND), allocatable, dimension(:,:) :: delsq_u
- real (kind=RKIND), allocatable, dimension(:,:) :: delsq_circulation, delsq_vorticity
-
-
- real (kind=RKIND), dimension(:,:), pointer :: u_src
- real (kind=RKIND), parameter :: rho_ref = 1000.0
-
- call timer_start("compute_tend_u")
-
- h => s % h % array
- u => s % u % array
- v => s % v % array
- wTop => s % wTop % array
- h_edge => s % h_edge % array
- circulation => s % circulation % array
- vorticity => s % vorticity % array
- divergence => s % divergence % array
- ke => s % ke % array
- ke_edge => s % ke_edge % array
- pv_edge => s % pv_edge % array
- MontPot => s % MontPot % array
- pressure => s % pressure % array
- vertViscTopOfEdge => d % vertViscTopOfEdge % array
-
- weightsOnEdge => grid % weightsOnEdge % array
- kiteAreasOnVertex => grid % kiteAreasOnVertex % array
- cellsOnEdge => grid % cellsOnEdge % array
- cellsOnVertex => grid % cellsOnVertex % array
- verticesOnEdge => grid % verticesOnEdge % array
- nEdgesOnCell => grid % nEdgesOnCell % array
- edgesOnCell => grid % edgesOnCell % array
- nEdgesOnEdge => grid % nEdgesOnEdge % array
- edgesOnEdge => grid % edgesOnEdge % array
- edgesOnVertex => grid % edgesOnVertex % array
- dcEdge => grid % dcEdge % array
- dvEdge => grid % dvEdge % array
- areaCell => grid % areaCell % array
- areaTriangle => grid % areaTriangle % array
- h_s => grid % h_s % array
-! mrp 110516 cleanup fvertex fedge not used in this subroutine
- fVertex => grid % fVertex % array
- fEdge => grid % fEdge % array
- zMidZLevel => grid % zMidZLevel % array
- zTopZLevel => grid % zTopZLevel % array
- maxLevelCell => grid % maxLevelCell % array
- maxLevelEdgeTop => grid % maxLevelEdgeTop % array
- maxLevelVertexBot => grid % maxLevelVertexBot % array
-
- tend_u => tend % u % array
-
- nCells = grid % nCells
- nEdges = grid % nEdges
- nEdgesSolve = grid % nEdgesSolve
- nVertices = grid % nVertices
- nVertLevels = grid % nVertLevels
-
- u_src => grid % u_src % array
-
- meshScalingDel2 => grid % meshScalingDel2 % array
- meshScalingDel4 => grid % meshScalingDel4 % array
-
- !
- ! velocity tendency: start accumulating tendency terms
- !
- ! mrp 110516 efficiency: could remove next line and have first tend_u operation not be additive
- tend_u(:,:) = 0.0
-
- !
- ! velocity tendency: nonlinear Coriolis term and grad of kinetic energy
- !
-
- call timer_start("compute_tend_u-coriolis")
-
-! call OcnVelCoriolisTend(grid, pv_edge, h_edge, u, ke, tend_u, err)
-
- do iEdge=1,grid % nEdgesSolve
- cell1 = cellsOnEdge(1,iEdge)
- cell2 = cellsOnEdge(2,iEdge)
-
- do k=1,maxLevelEdgeTop(iEdge)
-
- q = 0.0
- do j = 1,nEdgesOnEdge(iEdge)
- eoe = edgesOnEdge(j,iEdge)
- workpv = 0.5 * (pv_edge(k,iEdge) + pv_edge(k,eoe))
- q = q + weightsOnEdge(j,iEdge) * u(k,eoe) * workpv * h_edge(k,eoe)
- end do
-
- tend_u(k,iEdge) = tend_u(k,iEdge) &
- + q &
- - ( ke(k,cell2) - ke(k,cell1) ) / dcEdge(iEdge)
-
- end do
- end do
-
- call timer_stop("compute_tend_u-coriolis")
-
- !
- ! velocity tendency: vertical advection term -w du/dz
- !
- call timer_start("compute_tend_u-vert adv")
-
- call OcnVelVadvTend(grid, u, wTop, tend_u, err)
-
- call timer_stop("compute_tend_u-vert adv")
-
- !
- ! velocity tendency: pressure gradient
- !
- call timer_start("compute_tend_u-pressure grad")
-
- if (config_vert_grid_type.eq.'isopycnal') then
- call OcnVelPressureGradTend(grid, MontPot, tend_u, err)
- elseif (config_vert_grid_type.eq.'zlevel') then
- call OcnVelPressureGradTend(grid, pressure, tend_u, err)
- end if
-
- call timer_stop("compute_tend_u-pressure grad")
-
- !
- ! velocity tendency: del2 dissipation, </font>
<font color="black">u_2 </font>
<font color="red">abla^2 u
- ! computed as </font>
<font color="black">u( </font>
<font color="black">abla divergence + k \times </font>
<font color="gray">abla vorticity )
- ! strictly only valid for config_h_mom_eddy_visc2 == constant
- !
- call timer_start("compute_tend_u-horiz mix")
-
- call OcnVelHmixTend(grid, divergence, vorticity, tend_u, err)
-
- call timer_stop("compute_tend_u-horiz mix")
-
- !
- ! velocity tendency: forcing and bottom drag
- !
- ! mrp 101115 note: in order to include flux boundary conditions, we will need to
- ! know the bottom edge with nonzero velocity and place the drag there.
-
- call timer_start("compute_tend_u-forcings")
-
- call OcnVelForcingTend(grid, u, u_src, ke_edge, h_edge, tend_u, err)
-
- call timer_stop("compute_tend_u-forcings")
-
- !
- ! 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")
- call OcnVelVmixTendExplicit(grid, u, h_edge, vertViscTopOfEdge, tend_u, err)
-! 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
-
-! 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!}}}
-
subroutine filter_btr_mode_tend_u(tend, s, d, grid)!{{{
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Filter and remove barotropic mode from the tendencies
@@ -2158,887 +1685,6 @@
end subroutine filter_btr_mode_u!}}}
- subroutine compute_f_uperp(s, grid)!{{{
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
- ! Put f*uBcl^{perp} in u as a work variable
- !
- ! Input: s - current model state
- ! grid - grid metadata
- !
- ! Output: tend - computed tendencies for prognostic variables
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
- implicit none
-
- type (state_type), intent(inout) :: s
- type (mesh_type), intent(in) :: grid
-
-! mrp 110512 I just split compute_tend into compute_tend_u and compute_tend_h.
-! Some of these variables can be removed, but at a later time.
- integer :: iEdge, iCell, iVertex, k, cell1, cell2, &
- vertex1, vertex2, eoe, i, j
-
- integer :: nCells, nEdges, nVertices, nVertLevels, nEdgesSolve
- real (kind=RKIND) :: flux, vorticity_abs, h_vertex, workpv, q, &
- upstream_bias, wTopEdge, rho0Inv, r
- real (kind=RKIND), dimension(:), pointer :: &
- h_s, fVertex, fEdge, dvEdge, dcEdge, areaCell, areaTriangle, &
- zMidZLevel, zTopZLevel
- real (kind=RKIND), dimension(:,:), pointer :: &
- weightsOnEdge, kiteAreasOnVertex, h_edge, h, u, uBcl, v, pressure, &
- tend_u, circulation, vorticity, ke, ke_edge, pv_edge, &
- MontPot, wTop, divergence, vertViscTopOfEdge
- type (dm_info) :: dminfo
-
- integer, dimension(:), pointer :: nEdgesOnCell, nEdgesOnEdge, &
- maxLevelCell, maxLevelEdgeTop, maxLevelVertexBot
- integer, dimension(:,:), pointer :: &
- cellsOnEdge, cellsOnVertex, verticesOnEdge, edgesOnCell, &
- edgesOnEdge, edgesOnVertex
- real (kind=RKIND) :: u_diffusion
- real (kind=RKIND), dimension(:), allocatable:: fluxVertTop,w_dudzTopEdge
-
- real (kind=RKIND), allocatable, dimension(:,:) :: delsq_divergence
- real (kind=RKIND), allocatable, dimension(:,:) :: delsq_u
- real (kind=RKIND), allocatable, dimension(:,:) :: delsq_circulation, delsq_vorticity
-
-
- real (kind=RKIND), dimension(:,:), pointer :: u_src
- real (kind=RKIND), parameter :: rho_ref = 1000.0
-
- call timer_start("compute_f_uperp")
-
- h => s % h % array
- u => s % u % array
- uBcl => s % uBcl % array
- v => s % v % array
- wTop => s % wTop % array
- h_edge => s % h_edge % array
- circulation => s % circulation % array
- vorticity => s % vorticity % array
- divergence => s % divergence % array
- ke => s % ke % array
- ke_edge => s % ke_edge % array
- pv_edge => s % pv_edge % array
- MontPot => s % MontPot % array
- pressure => s % pressure % array
-
- weightsOnEdge => grid % weightsOnEdge % array
- kiteAreasOnVertex => grid % kiteAreasOnVertex % array
- cellsOnEdge => grid % cellsOnEdge % array
- cellsOnVertex => grid % cellsOnVertex % array
- verticesOnEdge => grid % verticesOnEdge % array
- nEdgesOnCell => grid % nEdgesOnCell % array
- edgesOnCell => grid % edgesOnCell % array
- nEdgesOnEdge => grid % nEdgesOnEdge % array
- edgesOnEdge => grid % edgesOnEdge % array
- edgesOnVertex => grid % edgesOnVertex % array
- dcEdge => grid % dcEdge % array
- dvEdge => grid % dvEdge % array
- areaCell => grid % areaCell % array
- areaTriangle => grid % areaTriangle % array
- h_s => grid % h_s % array
- fVertex => grid % fVertex % array
- fEdge => grid % fEdge % array
- zMidZLevel => grid % zMidZLevel % array
- zTopZLevel => grid % zTopZLevel % array
- maxLevelCell => grid % maxLevelCell % array
- maxLevelEdgeTop => grid % maxLevelEdgeTop % array
- maxLevelVertexBot => grid % maxLevelVertexBot % array
-
-
- nCells = grid % nCells
- nEdges = grid % nEdges
- nEdgesSolve = grid % nEdgesSolve
- nVertices = grid % nVertices
- nVertLevels = grid % nVertLevels
-
- !
- ! Put f*uBcl^{perp} in u as a work variable
- !
- do iEdge=1,grid % nEdgesSolve
- cell1 = cellsOnEdge(1,iEdge)
- cell2 = cellsOnEdge(2,iEdge)
-
- do k=1,maxLevelEdgeTop(iEdge)
-
- u(k,iEdge) = 0.0
- do j = 1,nEdgesOnEdge(iEdge)
- eoe = edgesOnEdge(j,iEdge)
- u(k,iEdge) = u(k,iEdge) + weightsOnEdge(j,iEdge) * uBcl(k,eoe) * fEdge(eoe)
- end do
- end do
- end do
-
- call timer_stop("compute_f_uperp")
-
- end subroutine compute_f_uperp!}}}
-
- subroutine compute_scalar_tend(tend, s, d, grid)!{{{
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
- !
- ! Input: s - current model state
- ! grid - grid metadata
- ! note: the variable s % tracers really contains the tracers,
- ! not tracers*h
- !
- ! Output: tend - computed scalar tendencies
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
- implicit none
-
- type (tend_type), intent(inout) :: tend
- type (state_type), intent(in) :: s
- type (diagnostics_type), intent(in) :: d
- type (mesh_type), intent(in) :: grid
-
- integer :: i, k, iCell, iEdge, iTracer, cell1, cell2, upwindCell,&
- nEdges, nCells, nCellsSolve, nVertLevels, num_tracers, err
- real (kind=RKIND) :: invAreaCell1, invAreaCell2, tracer_turb_flux
- real (kind=RKIND) :: flux, tracer_edge, r
- real (kind=RKIND), dimension(:), pointer :: &
- h_s, fVertex, fEdge, dvEdge, dcEdge, areaCell, areaTriangle
- real (kind=RKIND), dimension(:,:), pointer :: &
- u,h,wTop, h_edge, vertDiffTopOfCell
- real (kind=RKIND), dimension(:,:,:), pointer :: &
- tracers, tend_tr
- integer, dimension(:,:), pointer :: boundaryEdge
- type (dm_info) :: dminfo
-
- integer, dimension(:), pointer :: nEdgesOnCell, nEdgesOnEdge, &
- maxLevelCell, maxLevelEdgeTop, maxLevelVertexBot
- integer, dimension(:,:), pointer :: cellsOnEdge, boundaryCell
- real (kind=RKIND), dimension(:), pointer :: zTopZLevel,zMidZLevel, &
- hRatioZLevelK, hRatioZLevelKm1, meshScalingDel2, meshScalingDel4
- real (kind=RKIND), dimension(:), allocatable:: tracer2ndDer, tracersIn, tracersOut, posZMidZLevel, &
- posZTopZLevel
- real (kind=RKIND), dimension(:,:), allocatable:: fluxVertTop, boundaryMask
- real (kind=RKIND), dimension(:,:,:), allocatable::tr_flux, tr_div, delsq_tracer, tracerTop
-
-
- real (kind=RKIND) :: d2fdx2_cell1, d2fdx2_cell2
- real (kind=RKIND), dimension(:,:,:), pointer :: deriv_two
- real (kind=RKIND) :: coef_3rd_order, flux3Coef, cSignWTop
-
- integer :: index_temperature, index_salinity, rrr
- real (kind=RKIND), dimension(:), pointer :: temperatureRestore, salinityRestore
-
- call timer_start("compute_scalar_tend")
-
- u => s % u % array
- h => s % h % array
- boundaryCell=> grid % boundaryCell % array
- wTop => s % wTop % array
- tracers => s % tracers % array
- h_edge => s % h_edge % array
- vertDiffTopOfCell => d % vertDiffTopOfCell % array
-
- tend_tr => tend % tracers % array
-
- areaCell => grid % areaCell % array
- cellsOnEdge => grid % cellsOnEdge % array
- dvEdge => grid % dvEdge % array
- dcEdge => grid % dcEdge % array
- zTopZLevel => grid % zTopZLevel % array
- zMidZLevel => grid % zMidZLevel % array
- hRatioZLevelK => grid % hRatioZLevelK % array
- hRatioZLevelKm1 => grid % hRatioZLevelKm1 % array
- boundaryEdge => grid % boundaryEdge % array
- maxLevelCell => grid % maxLevelCell % array
- maxLevelEdgeTop => grid % maxLevelEdgeTop % array
- maxLevelVertexBot => grid % maxLevelVertexBot % array
-
- nEdges = grid % nEdges
- nCells = grid % nCells
- nCellsSolve = grid % nCellsSolve
- nVertLevels = grid % nVertLevels
- num_tracers = s % num_tracers
-
- meshScalingDel2 => grid % meshScalingDel2 % array
- meshScalingDel4 => grid % meshScalingDel4 % array
-
-
- deriv_two => grid % deriv_two % array
-
- !
- ! initialize tracer tendency (RHS of tracer equation) to zero.
- !
- tend_tr(:,:,:) = 0.0
-
- !
- ! tracer tendency: horizontal advection term -div( h \phi u)
- !
- ! mrp 101115 note: in order to include flux boundary conditions, we will need to
- ! assign h_edge for maxLevelEdgeTop:maxLevelEdgeBot in the compute_solve_diagnostics
- ! and then change maxLevelEdgeTop to maxLevelEdgeBot in the following section.
- ! tracer_edge at the boundary will also need to be defined for flux boundaries.
-
- call timer_start("compute_scalar_tend-horiz adv")
-
- call OcnTracerHadvTend(grid, u, h_edge, tracers, tend_tr, err)
-
- call timer_stop("compute_scalar_tend-horiz adv")
-
-
- !
- ! tracer tendency: vertical advection term -d/dz( h \phi w)
- !
-
- call timer_start("compute_scalar_tend-vert adv")
-
- call OcnTracerVadvTend(grid, wTop, tracers, tend_tr, err)
-
- call timer_stop("compute_scalar_tend-vert adv")
-
- !
- ! tracer tendency: del2 horizontal tracer diffusion, div(h \kappa_2 </font>
<font color="red">abla \phi)
- !
- call timer_start("compute_scalar_tend-horiz diff")
-
- call OcnTracerHmixTend(grid, h_edge, tracers, tend_tr, err)
-
- call timer_stop("compute_scalar_tend-horiz diff")
-
-! mrp 110516 printing
-!print *, 'tend_tr 1',minval(tend_tr(3,1,1:nCells)),&
-! maxval(tend_tr(3,1,1:nCells))
-!print *, 'tracer 1',minval(tracers(3,1,1:nCells)),&
-! maxval(tracers(3,1,1:nCells))
-! mrp 110516 printing end
-
- !
- ! 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")
- call OcnTracerVmixTendExplicit(grid, h, vertDiffTopOfCell, tracers, tend_tr, err)
-! 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_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))
-! mrp 110516 printing end
-
- !
- ! add restoring to T and S in top model layer
- !
- call timer_start("compute_scalar_tend-restoring")
-
- call OcnRestoringTend(grid, h, s%index_temperature, s%index_salinity, tracers, tend_tr, err)
-
- call timer_stop("compute_scalar_tend-restoring")
-
- 10 format(2i8,10e20.10)
- call timer_stop("compute_scalar_tend")
-
- end subroutine compute_scalar_tend!}}}
-
- subroutine compute_solve_diagnostics(dt, s, grid)!{{{
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
- ! Compute diagnostic fields used in the tendency computations
- !
- ! Input: grid - grid metadata
- !
- ! Output: s - computed diagnostics
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
- implicit none
-
- real (kind=RKIND), intent(in) :: dt
- type (state_type), intent(inout) :: s
- type (mesh_type), intent(in) :: grid
-
-
- integer :: iEdge, iCell, iVertex, k, cell1, cell2, vertex1, vertex2, eoe, i, j, cov
- real (kind=RKIND) :: flux, vorticity_abs, h_vertex, workpv, rho0Inv
-
- integer :: nCells, nEdges, nVertices, nVertLevels, vertexDegree, fCoef
-
-
- real (kind=RKIND), dimension(:), pointer :: &
- h_s, fVertex, fEdge, dvEdge, dcEdge, areaCell, areaTriangle, &
- hZLevel
- real (kind=RKIND), dimension(:,:), pointer :: &
- weightsOnEdge, kiteAreasOnVertex, h_edge, h, u, v, w, pressure,&
- circulation, vorticity, ke, ke_edge, MontPot, wTop, &
- pv_edge, pv_vertex, pv_cell, gradPVn, gradPVt, divergence, &
- rho, temperature, salinity
- real (kind=RKIND), dimension(:,:,:), pointer :: tracers
- real (kind=RKIND), dimension(:), allocatable:: pTop
- real (kind=RKIND), dimension(:,:), allocatable:: div_u
- character :: c1*6
-
- integer, dimension(:,:), pointer :: cellsOnEdge, cellsOnVertex, &
- verticesOnEdge, edgesOnCell, edgesOnEdge, edgesOnVertex, &
- boundaryEdge, boundaryCell
- integer, dimension(:), pointer :: nEdgesOnCell, nEdgesOnEdge, &
- maxLevelCell, maxLevelEdgeTop, maxLevelEdgeBot, &
- maxLevelVertexBot, maxLevelVertexTop
- real (kind=RKIND) :: d2fdx2_cell1, d2fdx2_cell2
- real (kind=RKIND), dimension(:,:,:), pointer :: deriv_two
- real (kind=RKIND) :: coef_3rd_order
- real (kind=RKIND) :: r, h1, h2
-
- call timer_start("compute_solve_diagnostics")
-
- h => s % h % array
- u => s % u % array
- v => s % v % array
- wTop => s % wTop % array
- h_edge => s % h_edge % array
- circulation => s % circulation % array
- vorticity => s % vorticity % array
- divergence => s % divergence % array
- ke => s % ke % array
- ke_edge => s % ke_edge % array
- pv_edge => s % pv_edge % array
- pv_vertex => s % pv_vertex % array
- pv_cell => s % pv_cell % array
- gradPVn => s % gradPVn % array
- gradPVt => s % gradPVt % array
- rho => s % rho % array
- tracers => s % tracers % array
- MontPot => s % MontPot % array
- pressure => s % pressure % array
-
- weightsOnEdge => grid % weightsOnEdge % array
- kiteAreasOnVertex => grid % kiteAreasOnVertex % array
- cellsOnEdge => grid % cellsOnEdge % array
- cellsOnVertex => grid % cellsOnVertex % array
- verticesOnEdge => grid % verticesOnEdge % array
- nEdgesOnCell => grid % nEdgesOnCell % array
- edgesOnCell => grid % edgesOnCell % array
- nEdgesOnEdge => grid % nEdgesOnEdge % array
- edgesOnEdge => grid % edgesOnEdge % array
- edgesOnVertex => grid % edgesOnVertex % array
- dcEdge => grid % dcEdge % array
- dvEdge => grid % dvEdge % array
- areaCell => grid % areaCell % array
- areaTriangle => grid % areaTriangle % array
- h_s => grid % h_s % array
- fVertex => grid % fVertex % array
- fEdge => grid % fEdge % array
- hZLevel => grid % hZLevel % array
- deriv_two => grid % deriv_two % array
- maxLevelCell => grid % maxLevelCell % array
- maxLevelEdgeTop => grid % maxLevelEdgeTop % array
- maxLevelEdgeBot => grid % maxLevelEdgeBot % array
- maxLevelVertexBot => grid % maxLevelVertexBot % array
- maxLevelVertexTop => grid % maxLevelVertexTop % array
-
- nCells = grid % nCells
- nEdges = grid % nEdges
- nVertices = grid % nVertices
- nVertLevels = grid % nVertLevels
- vertexDegree = grid % vertexDegree
-
- boundaryEdge => grid % boundaryEdge % array
- boundaryCell => grid % boundaryCell % array
-
- !
- ! Compute height on cell edges at velocity locations
- ! Namelist options control the order of accuracy of the reconstructed h_edge value
- !
- ! mrp 101115 note: in order to include flux boundary conditions, we will need to
- ! assign h_edge for maxLevelEdgeTop:maxLevelEdgeBot in the following section
-
- ! mrp 110516 efficiency note: For z-level, only do this on level 1. h_edge for all
- ! lower levels is defined by hZlevel.
-
- call timer_start("compute_solve_diagnostics-hEdge")
-
- coef_3rd_order = 0.
- if (config_thickness_adv_order == 3) coef_3rd_order = 1.0
- if (config_thickness_adv_order == 3 .and. config_monotonic) coef_3rd_order = 0.25
-
- if (config_thickness_adv_order == 2) then
- call timer_start("compute_solve_diagnostics-hEdge 2")
-
- do iEdge=1,nEdges
- 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
- end do
- call timer_stop("compute_solve_diagnostics-hEdge 2")
-
- else if (config_thickness_adv_order == 3) then
- call timer_start("compute_solve_diagnostics-hEdge 3")
-
- do iEdge=1,nEdges
- cell1 = cellsOnEdge(1,iEdge)
- cell2 = cellsOnEdge(2,iEdge)
-
- do k=1,maxLevelEdgeTop(iEdge)
-
- d2fdx2_cell1 = 0.0
- d2fdx2_cell2 = 0.0
-
- !-- if not a boundary cell
- if(boundaryCell(k,cell1).eq.0.and.boundaryCell(k,cell2).eq.0) then
-
- d2fdx2_cell1 = deriv_two(1,1,iEdge) * h(k,cell1)
- d2fdx2_cell2 = deriv_two(1,2,iEdge) * h(k,cell2)
-
- !-- all edges of cell 1
- do i=1, grid % nEdgesOnCell % array (cell1)
- d2fdx2_cell1 = d2fdx2_cell1 + &
- deriv_two(i+1,1,iEdge) * h(k,grid % CellsOnCell % array (i,cell1))
- end do
-
- !-- all edges of cell 2
- do i=1, grid % nEdgesOnCell % array (cell2)
- d2fdx2_cell2 = d2fdx2_cell2 + &
- deriv_two(i+1,2,iEdge) * h(k,grid % CellsOnCell % array (i,cell2))
- end do
-
- endif
-
- !-- if u > 0:
- if (u(k,iEdge) > 0) then
- h_edge(k,iEdge) = &
- 0.5*(h(k,cell1) + h(k,cell2)) &
- -(dcEdge(iEdge) **2) * (d2fdx2_cell1 + d2fdx2_cell2) / 12. &
- -(dcEdge(iEdge) **2) * coef_3rd_order*(d2fdx2_cell1 - d2fdx2_cell2) / 12.
- !-- else u <= 0:
- else
- h_edge(k,iEdge) = &
- 0.5*(h(k,cell1) + h(k,cell2)) &
- -(dcEdge(iEdge) **2) * (d2fdx2_cell1 + d2fdx2_cell2) / 12. &
- +(dcEdge(iEdge) **2) * coef_3rd_order*(d2fdx2_cell1 - d2fdx2_cell2) / 12.
- end if
-
- end do ! do k
- end do ! do iEdge
-
- call timer_stop("compute_solve_diagnostics-hEdge 3")
- else if (config_thickness_adv_order == 4) then
- call timer_start("compute_solve_diagnostics-hEdge 4")
-
- do iEdge=1,nEdges
- cell1 = cellsOnEdge(1,iEdge)
- cell2 = cellsOnEdge(2,iEdge)
-
- do k=1,maxLevelEdgeTop(iEdge)
-
- d2fdx2_cell1 = 0.0
- d2fdx2_cell2 = 0.0
-
- !-- if not a boundary cell
- if(boundaryCell(k,cell1).eq.0.and.boundaryCell(k,cell2).eq.0) then
-
- d2fdx2_cell1 = deriv_two(1,1,iEdge) * h(k,cell1)
- d2fdx2_cell2 = deriv_two(1,2,iEdge) * h(k,cell2)
-
- !-- all edges of cell 1
- do i=1, grid % nEdgesOnCell % array (cell1)
- d2fdx2_cell1 = d2fdx2_cell1 + &
- deriv_two(i+1,1,iEdge) * h(k,grid % CellsOnCell % array (i,cell1))
- end do
-
- !-- all edges of cell 2
- do i=1, grid % nEdgesOnCell % array (cell2)
- d2fdx2_cell2 = d2fdx2_cell2 + &
- deriv_two(i+1,2,iEdge) * h(k,grid % CellsOnCell % array (i,cell2))
- end do
-
- endif
-
- h_edge(k,iEdge) = &
- 0.5*(h(k,cell1) + h(k,cell2)) &
- -(dcEdge(iEdge) **2) * (d2fdx2_cell1 + d2fdx2_cell2) / 12.
-
- end do ! do k
- end do ! do iEdge
-
- call timer_stop("compute_solve_diagnostics-hEdge 4")
- endif ! if(config_thickness_adv_order == 2)
- call timer_stop("compute_solve_diagnostics-hEdge")
-
- !
- ! set the velocity and height at dummy address
- ! used -1e34 so error clearly occurs if these values are used.
- !
-!mrp 110516 change to zero, change back later:
- u(:,nEdges+1) = -1e34
- h(:,nCells+1) = -1e34
- tracers(s % index_temperature,:,nCells+1) = -1e34
- tracers(s % index_salinity,:,nCells+1) = -1e34
-
- !
- ! Compute circulation and relative vorticity at each vertex
- !
- circulation(:,:) = 0.0
- do iEdge=1,nEdges
- vertex1 = verticesOnEdge(1,iEdge)
- vertex2 = verticesOnEdge(2,iEdge)
- do k=1,maxLevelEdgeBot(iEdge)
- circulation(k,vertex1) = circulation(k,vertex1) - dcEdge(iEdge) * u(k,iEdge)
- circulation(k,vertex2) = circulation(k,vertex2) + dcEdge(iEdge) * u(k,iEdge)
- end do
- end do
- do iVertex=1,nVertices
- do k=1,maxLevelVertexBot(iVertex)
- vorticity(k,iVertex) = circulation(k,iVertex) / areaTriangle(iVertex)
- end do
- end do
-
- !
- ! Compute the divergence at each cell center
- !
- divergence(:,:) = 0.0
- do iEdge=1,nEdges
- cell1 = cellsOnEdge(1,iEdge)
- cell2 = cellsOnEdge(2,iEdge)
- do k=1,maxLevelEdgeBot(iEdge)
- divergence(k,cell1) = divergence(k,cell1) + u(k,iEdge)*dvEdge(iEdge)
- divergence(k,cell2) = divergence(k,cell2) - u(k,iEdge)*dvEdge(iEdge)
- enddo
- end do
- do iCell = 1,nCells
- r = 1.0 / areaCell(iCell)
- do k = 1,maxLevelCell(iCell)
- divergence(k,iCell) = divergence(k,iCell) * r
- enddo
- enddo
-
- !
- ! Compute kinetic energy in each cell
- !
- ke(:,:) = 0.0
- do iEdge=1,nEdges
- cell1 = cellsOnEdge(1,iEdge)
- cell2 = cellsOnEdge(2,iEdge)
- do k=1,maxLevelEdgeBot(iEdge)
- ke(k,cell1) = ke(k,cell1) + 0.25 * dcEdge(iEdge) * dvEdge(iEdge) * u(k,iEdge)**2.0
- ke(k,cell2) = ke(k,cell2) + 0.25 * dcEdge(iEdge) * dvEdge(iEdge) * u(k,iEdge)**2.0
- enddo
- end do
- do iCell = 1,nCells
- do k = 1,maxLevelCell(iCell)
- ke(k,iCell) = ke(k,iCell) / areaCell(iCell)
- enddo
- enddo
-
- !
- ! Compute v (tangential) velocities
- !
- v(:,:) = 0.0
- do iEdge = 1,nEdges
- do i=1,nEdgesOnEdge(iEdge)
- eoe = edgesOnEdge(i,iEdge)
- ! mrp 101115 note: in order to include flux boundary conditions,
- ! the following loop may need to change to maxLevelEdgeBot
- do k = 1,maxLevelEdgeTop(iEdge)
- v(k,iEdge) = v(k,iEdge) + weightsOnEdge(i,iEdge) * u(k, eoe)
- end do
- end do
- end do
-
- !
- ! Compute ke on cell edges at velocity locations for quadratic bottom drag.
- !
- ! mrp 101025 efficiency note: we could get rid of ke_edge completely by
- ! using sqrt(u(k,iEdge)**2 + v(k,iEdge)**2) in its place elsewhere.
- ke_edge = 0.0 !mrp remove 0 for efficiency
- do iEdge=1,nEdges
- cell1 = cellsOnEdge(1,iEdge)
- cell2 = cellsOnEdge(2,iEdge)
- do k=1,maxLevelEdgeTop(iEdge)
- ke_edge(k,iEdge) = 0.5 * (ke(k,cell1) + ke(k,cell2))
- end do
- end do
-
- !
- ! Compute height at vertices, pv at vertices, and average pv to edge locations
- ! ( this computes pv_vertex at all vertices bounding real cells and distance-1 ghost cells )
- !
- if (trim(config_time_integration) == 'RK4') then
- ! for RK4, PV is really PV = (eta+f)/h
- fCoef = 1
- elseif (trim(config_time_integration) == 'split_explicit' &
- .or.trim(config_time_integration) == 'unsplit_explicit') then
- ! for split explicit, PV is eta/h because f is added separately to the momentum forcing.
-! mrp temp, new should be:
- fCoef = 0
-! old, for testing:
-! fCoef = 1
- end if
-
- do iVertex = 1,nVertices
- do k=1,maxLevelVertexBot(iVertex)
- h_vertex = 0.0
- do i=1,vertexDegree
- h_vertex = h_vertex + h(k,cellsOnVertex(i,iVertex)) * kiteAreasOnVertex(i,iVertex)
- end do
- h_vertex = h_vertex / areaTriangle(iVertex)
-
- pv_vertex(k,iVertex) = (fCoef*fVertex(iVertex) + vorticity(k,iVertex)) / h_vertex
- end do
- end do
-
- !
- ! Compute pv at cell centers
- ! ( this computes pv_cell for all real cells and distance-1 ghost cells )
- !
- pv_cell(:,:) = 0.0
- do iVertex = 1,nVertices
- do i=1,vertexDegree
- iCell = cellsOnVertex(i,iVertex)
- do k = 1,maxLevelCell(iCell)
- pv_cell(k,iCell) = pv_cell(k,iCell) &
- + kiteAreasOnVertex(i, iVertex) * pv_vertex(k, iVertex) &
- / areaCell(iCell)
- enddo
- enddo
- enddo
-
- !
- ! Compute pv at the edges
- ! ( this computes pv_edge at all edges bounding real cells )
- !
- pv_edge(:,:) = 0.0
- do iVertex = 1,nVertices
- do i=1,vertexDegree
- iEdge = edgesOnVertex(i,iVertex)
- do k=1,maxLevelEdgeBot(iEdge)
- pv_edge(k,iEdge) = pv_edge(k,iEdge) + 0.5 * pv_vertex(k,iVertex)
- enddo
- end do
- end do
-
- !
- ! Compute gradient of PV in normal direction
- ! ( this computes gradPVn for all edges bounding real cells )
- !
- gradPVn(:,:) = 0.0
- do iEdge = 1,nEdges
- do k=1,maxLevelEdgeTop(iEdge)
- gradPVn(k,iEdge) = ( pv_cell(k,cellsOnEdge(2,iEdge)) &
- - pv_cell(k,cellsOnEdge(1,iEdge))) &
- / dcEdge(iEdge)
- enddo
- enddo
-
- !
- ! Compute gradient of PV in the tangent direction
- ! ( this computes gradPVt at all edges bounding real cells and distance-1 ghost cells )
- !
- do iEdge = 1,nEdges
- do k = 1,maxLevelEdgeBot(iEdge)
- gradPVt(k,iEdge) = ( pv_vertex(k,verticesOnEdge(2,iEdge)) &
- - pv_vertex(k,verticesOnEdge(1,iEdge))) &
- /dvEdge(iEdge)
- enddo
- enddo
-
- !
- ! Modify PV edge with upstream bias.
- !
- do iEdge = 1,nEdges
- do k = 1,maxLevelEdgeBot(iEdge)
- pv_edge(k,iEdge) = pv_edge(k,iEdge) &
- - 0.5 * dt* ( u(k,iEdge) * gradPVn(k,iEdge) &
- + v(k,iEdge) * gradPVt(k,iEdge) )
- enddo
- enddo
-
- !
- ! equation of state
- !
- ! For an isopycnal model, density should remain constant.
- ! For zlevel, calculate in-situ density
- 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
- !
- if (config_vert_grid_type.eq.'isopycnal') then
-
- ! For Isopycnal model.
- ! Compute pressure at top of each layer, and then
- ! Montgomery Potential.
- allocate(pTop(nVertLevels))
- do iCell=1,nCells
-
- ! assume atmospheric pressure at the surface is zero for now.
- pTop(1) = 0.0
- ! For isopycnal mode, p is the Montgomery Potential.
- ! At top layer it is g*SSH, where SSH may be off by a
- ! constant (ie, h_s can be relative to top or bottom)
- MontPot(1,iCell) = gravity &
- * (h_s(iCell) + sum(h(1:nVertLevels,iCell)))
-
- do k=2,nVertLevels
- pTop(k) = pTop(k-1) + rho(k-1,iCell)*gravity* h(k-1,iCell)
-
- ! from delta M = p delta / rho
- MontPot(k,iCell) = MontPot(k-1,iCell) &
- + pTop(k)*(1.0/rho(k,iCell) - 1.0/rho(k-1,iCell))
- end do
-
- end do
- deallocate(pTop)
-
- elseif (config_vert_grid_type.eq.'zlevel') then
-
- ! For z-level model.
- ! Compute pressure at middle of each level.
- ! At k=1, where p is pressure at a depth of hZLevel(1)/2, not
- ! pressure at middle of layer including SSH.
-
- do iCell=1,nCells
- ! compute pressure for z-level coordinates
- ! assume atmospheric pressure at the surface is zero for now.
-
- pressure(1,iCell) = rho(1,iCell)*gravity &
- * (h(1,iCell)-0.5*hZLevel(1))
-
- do k=2,maxLevelCell(iCell)
- pressure(k,iCell) = pressure(k-1,iCell) &
- + 0.5*gravity*( rho(k-1,iCell)*hZLevel(k-1) &
- + rho(k ,iCell)*hZLevel(k ))
- end do
-
- end do
-
- endif
-
- call compute_wtop(s,grid)
-
- call timer_stop("compute_solve_diagnostics")
-
- end subroutine compute_solve_diagnostics!}}}
-
- subroutine compute_wtop(s, grid)!{{{
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
- ! Compute diagnostic fields used in the tendency computations
- !
- ! Input: grid - grid metadata
- !
- ! Output: s - computed diagnostics
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
- implicit none
-
- type (state_type), intent(inout) :: s
- type (mesh_type), intent(in) :: grid
-
- ! mrp 110512 could clean this out, remove pointers?
- integer :: iEdge, iCell, iVertex, k, cell1, cell2, vertex1, vertex2, eoe, i, j, cov
- real (kind=RKIND) :: flux, vorticity_abs, h_vertex, workpv, rho0Inv
-
- integer :: nCells, nEdges, nVertices, nVertLevels, vertexDegree
-
-
- real (kind=RKIND), dimension(:), pointer :: &
- h_s, fVertex, fEdge, dvEdge, dcEdge, areaCell, areaTriangle, &
- hZLevel
- real (kind=RKIND), dimension(:,:), pointer :: u,wTop
- real (kind=RKIND), dimension(:,:), allocatable:: div_u
-
- integer, dimension(:,:), pointer :: cellsOnEdge, cellsOnVertex, &
- verticesOnEdge, edgesOnCell, edgesOnEdge, edgesOnVertex, &
- boundaryEdge, boundaryCell
- integer, dimension(:), pointer :: nEdgesOnCell, nEdgesOnEdge, &
- maxLevelCell, maxLevelEdgeTop, maxLevelEdgeBot, &
- maxLevelVertexBot, maxLevelVertexTop
-
- call timer_start("wTop")
-
- u => s % u % array
- wTop => s % wTop % array
-
- areaCell => grid % areaCell % array
- cellsOnEdge => grid % cellsOnEdge % array
- hZLevel => grid % hZLevel % array
- maxLevelCell => grid % maxLevelCell % array
- maxLevelEdgeBot => grid % maxLevelEdgeBot % array
- dvEdge => grid % dvEdge % array
-
- nCells = grid % nCells
- nEdges = grid % nEdges
- nVertLevels = grid % nVertLevels
-
- !
- ! vertical velocity through layer interface
- !
- if (config_vert_grid_type.eq.'isopycnal') then
- ! set vertical velocity to zero in isopycnal case
- wTop=0.0
-
- elseif (config_vert_grid_type.eq.'zlevel') then
-
- !
- ! Compute div(u) for each cell
- ! See Ringler et al. (2010) jcp paper, eqn 19, 21, and fig. 3.
- !
- allocate(div_u(nVertLevels,nCells+1))
- div_u(:,:) = 0.0
- do iEdge=1,nEdges
- cell1 = cellsOnEdge(1,iEdge)
- cell2 = cellsOnEdge(2,iEdge)
- do k=2,maxLevelEdgeBot(iEdge)
- flux = u(k,iEdge) * dvEdge(iEdge)
- div_u(k,cell1) = div_u(k,cell1) + flux
- div_u(k,cell2) = div_u(k,cell2) - flux
- end do
- end do
-
- do iCell=1,nCells
- ! Vertical velocity through layer interface at top and
- ! bottom is zero.
- wTop(1,iCell) = 0.0
- wTop(maxLevelCell(iCell)+1,iCell) = 0.0
- do k=maxLevelCell(iCell),2,-1
- wTop(k,iCell) = wTop(k+1,iCell) &
- - div_u(k,iCell)/areaCell(iCell)*hZLevel(k)
- end do
- end do
- deallocate(div_u)
-
- endif
-
- call timer_stop("wTop")
-
- end subroutine compute_wtop!}}}
-
subroutine enforce_boundaryEdge(tend, grid)!{{{
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Enforce any boundary conditions on the normal velocity at each edge
</font>
</pre>