<p><b>dwj07@fsu.edu</b> 2011-08-31 11:56:45 -0600 (Wed, 31 Aug 2011)</p><p><br>
Modularizing horizontal advection of tracers.<br>
<br>
Fixing some comments.<br>
<br>
Fixing a bug that in HmixTracerDel2.<br>
</p><hr noshade><pre><font color="gray">Modified: branches/ocean_projects/performance/mpas/src/core_ocean/Makefile
===================================================================
--- branches/ocean_projects/performance/mpas/src/core_ocean/Makefile 2011-08-30 22:18:11 UTC (rev 968)
+++ branches/ocean_projects/performance/mpas/src/core_ocean/Makefile 2011-08-31 17:56:45 UTC (rev 969)
@@ -15,6 +15,10 @@
module_OcnVmixVel.o \
module_OcnVmixVelConst.o \
module_OcnVmixVelTanh.o \
+ module_OcnHorizAdvTracer.o \
+ module_OcnHorizAdvTracer2nd.o \
+ module_OcnHorizAdvTracer3rd.o \
+ module_OcnHorizAdvTracer4th.o \
module_OcnCoriolis.o \
module_OcnPressureGrad.o \
module_OcnThickness.o \
@@ -36,6 +40,14 @@
module_advection.o:
+module_OcnHorizAdvTracer2nd.o :
+
+module_OcnHorizAdvTracer3rd.o :
+
+module_OcnHorizAdvTracer4th.o :
+
+module_OcnHorizAdvTracer.o : module_OcnHorizAdvTracer2nd.o module_OcnHorizAdvTracer3rd.o module_OcnHorizAdvTracer4th.o
+
module_OcnPressureGrad.o:
module_OcnCoriolis.o:
@@ -70,7 +82,31 @@
module_global_diagnostics.o:
-module_mpas_core.o: module_advection.o module_global_diagnostics.o module_test_cases.o module_OcnVelVertAdv.o module_OcnPressureGrad.o module_OcnVelocityForcing.o module_OcnThickness.o module_OcnCoriolis.o module_OcnVmixVel.o module_OcnVmixVelConst.o module_OcnVmixVelTanh.o module_OcnVmixTracer.o module_OcnVmixTracerConst.o module_OcnVmixTracerTanh.o module_OcnHmixTracer.o module_OcnHmixTracerDel2.o module_OcnHmixTracerDel4.o module_OcnHmixVel.o module_OcnHmixVelDel2.o module_OcnHmixVelDel4.o module_time_integration.o
+module_mpas_core.o: module_advection.o \
+ module_global_diagnostics.o \
+ module_test_cases.o \
+ module_OcnHorizAdvTracer2nd.o \
+ module_OcnHorizAdvTracer3rd.o \
+ module_OcnHorizAdvTracer4th.o \
+ module_OcnHorizAdvTracer.o \
+ module_OcnVelVertAdv.o \
+ module_OcnPressureGrad.o \
+ module_OcnVelocityForcing.o \
+ module_OcnThickness.o \
+ module_OcnCoriolis.o \
+ module_OcnVmixVel.o \
+ module_OcnVmixVelConst.o \
+ module_OcnVmixVelTanh.o \
+ module_OcnVmixTracer.o \
+ module_OcnVmixTracerConst.o \
+ module_OcnVmixTracerTanh.o \
+ module_OcnHmixTracer.o \
+ module_OcnHmixTracerDel2.o \
+ module_OcnHmixTracerDel4.o \
+ module_OcnHmixVel.o \
+ module_OcnHmixVelDel2.o \
+ module_OcnHmixVelDel4.o \
+ module_time_integration.o
clean:
$(RM) *.o *.mod *.f90 libdycore.a
Modified: branches/ocean_projects/performance/mpas/src/core_ocean/module_OcnHmixTracerDel2.F
===================================================================
--- branches/ocean_projects/performance/mpas/src/core_ocean/module_OcnHmixTracerDel2.F 2011-08-30 22:18:11 UTC (rev 968)
+++ branches/ocean_projects/performance/mpas/src/core_ocean/module_OcnHmixTracerDel2.F 2011-08-31 17:56:45 UTC (rev 969)
@@ -170,7 +170,7 @@
!*** compute a boundary mask at each vertical level
!*** on this edge
- boundaryMask = boundaryEdge(k,iEdge) * 1.d0
+ boundaryMask = (.not.boundaryEdge(k,iEdge)) * 1.d0
do iTracer=1,nTracers
Modified: branches/ocean_projects/performance/mpas/src/core_ocean/module_OcnHmixTracerDel4.F
===================================================================
--- branches/ocean_projects/performance/mpas/src/core_ocean/module_OcnHmixTracerDel4.F 2011-08-30 22:18:11 UTC (rev 968)
+++ branches/ocean_projects/performance/mpas/src/core_ocean/module_OcnHmixTracerDel4.F 2011-08-31 17:56:45 UTC (rev 969)
@@ -174,7 +174,7 @@
!*** compute a boundary mask at each vertical level
!*** on this edge
- boundaryMask = grid % boundaryEdge % array(k,iEdge) * 1.d0
+ boundaryMask = (.not.grid % boundaryEdge % array(k,iEdge)) * 1.d0
do iTracer=1,nTracers
Added: branches/ocean_projects/performance/mpas/src/core_ocean/module_OcnHorizAdvTracer.F
===================================================================
--- branches/ocean_projects/performance/mpas/src/core_ocean/module_OcnHorizAdvTracer.F (rev 0)
+++ branches/ocean_projects/performance/mpas/src/core_ocean/module_OcnHorizAdvTracer.F 2011-08-31 17:56:45 UTC (rev 969)
@@ -0,0 +1,166 @@
+!|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
+!
+! OcnHorizAdvTracer
+!
+!> \brief MPAS ocean horizontal tracer mixing driver
+!> \author Doug Jacobsen
+!> \date 31 August 2011
+!> \version SVN:$Id:$
+!> \details
+!> This module contains the main driver routine for computing
+!> horizontal advection tendencies for tracers.
+!>
+!>
+!
+!-----------------------------------------------------------------------
+
+module OcnHorizAdvTracer
+
+ use grid_types
+ use configure
+ use OcnHorizAdvTracer2nd
+ use OcnHorizAdvTracer3rd
+ use OcnHorizAdvTracer4th
+
+ implicit none
+ private
+ save
+
+ !--------------------------------------------------------------------
+ !
+ ! Public parameters
+ !
+ !--------------------------------------------------------------------
+
+ !--------------------------------------------------------------------
+ !
+ ! Public member functions
+ !
+ !--------------------------------------------------------------------
+
+ public :: OcnHorizAdvTracerTend, &
+ OcnHorizAdvTracerInit
+
+ !--------------------------------------------------------------------
+ !
+ ! Private module variables
+ !
+ !--------------------------------------------------------------------
+
+
+!***********************************************************************
+
+contains
+
+!***********************************************************************
+!
+! routine OcnHorizAdvTracerTend
+!
+!> \brief Computes tendency term for horizontal tracer mixing
+!> \author Doug Jacobsen
+!> \date 31 August 2011
+!> \version SVN:$Id$
+!> \details
+!> This routine computes the horizontal advection tendency for tracers
+!> based on current state and user choices of advection order.
+!>
+!
+!-----------------------------------------------------------------------
+
+ subroutine OcnHorizAdvTracerTend(tend, tracers, hEdge, u, grid)
+
+ !-----------------------------------------------------------------
+ !
+ ! input variables
+ !
+ !-----------------------------------------------------------------
+
+ real (kind=RKIND), dimension(:,:,:), intent(in) :: &
+ tracers !< Input: array of tracers
+
+ real (kind=RKIND), dimension(:,:), intent(in) :: &
+ hEdge !< Input: layer thickness at cell edges
+
+ real (kind=RKIND), dimension(:,:), intent(in) :: &
+ u !< Input: Velocity
+
+ type (mesh_type), intent(in) :: &
+ grid !< Input: grid information
+
+ !-----------------------------------------------------------------
+ !
+ ! input/output variables
+ !
+ !-----------------------------------------------------------------
+
+ real (kind=RKIND), dimension(:,:,:), intent(inout) :: &
+ tend !< Input/Output: tracer tendency
+
+ !-----------------------------------------------------------------
+ !
+ ! output variables
+ !
+ !-----------------------------------------------------------------
+
+ !-----------------------------------------------------------------
+ !
+ ! local variables
+ !
+ !-----------------------------------------------------------------
+
+ !-----------------------------------------------------------------
+ !
+ ! call relevant routines for computing tendencies
+ ! note that the user can choose multiple options and the
+ ! tendencies will be added together
+ !
+ !-----------------------------------------------------------------
+
+ call OcnHorizAdvTracer2ndTend(tend, tracers, hEdge, u, grid)
+ call OcnHorizAdvTracer3rdTend(tend, tracers, hEdge, u, grid)
+ call OcnHorizAdvTracer4thTend(tend, tracers, hEdge, u, grid)
+
+ !--------------------------------------------------------------------
+
+ end subroutine OcnHorizAdvTracerTend
+
+!***********************************************************************
+!
+! routine OcnHorizAdvTracerInit
+!
+!> \brief Initializes ocean tracer horizontal mixing quantities
+!> \author Doug Jacobsen
+!> \date 31 August 2011
+!> \version SVN:$Id$
+!> \details
+!> This routine initializes a variety of quantities related to
+!> horizontal tracer advection in the ocean. Since a variety of
+!> orders of accuracy are available, this routine primarily calls the
+!> individual init routines for each parameterization.
+!
+!-----------------------------------------------------------------------
+
+
+ subroutine OcnHorizAdvTracerInit
+
+ !--------------------------------------------------------------------
+
+ !-----------------------------------------------------------------
+ !
+ ! call individual init routines for each parameterization
+ !
+ !-----------------------------------------------------------------
+
+ call OcnHorizAdvTracer2ndInit
+ call OcnHorizAdvTracer3rdInit
+ call OcnHorizAdvTracer4thInit
+
+ !--------------------------------------------------------------------
+
+ end subroutine OcnHorizAdvTracerInit
+
+!***********************************************************************
+
+end module OcnHorizAdvTracer
+
+!|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
Added: branches/ocean_projects/performance/mpas/src/core_ocean/module_OcnHorizAdvTracer2nd.F
===================================================================
--- branches/ocean_projects/performance/mpas/src/core_ocean/module_OcnHorizAdvTracer2nd.F (rev 0)
+++ branches/ocean_projects/performance/mpas/src/core_ocean/module_OcnHorizAdvTracer2nd.F 2011-08-31 17:56:45 UTC (rev 969)
@@ -0,0 +1,201 @@
+!|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
+!
+! OcnHorizAdvTracer2nd
+!
+!> \brief Ocean horizontal advection - 2nd order
+!> \author Doug Jacobsen
+!> \date 31 August 2011
+!> \version SVN:$Id:$
+!> \details
+!> This module contains routines for computing horizontal advection
+!> tendencies using a 2nd order formulation.
+!
+!-----------------------------------------------------------------------
+
+module OcnHorizAdvTracer2nd
+
+ use grid_types
+ use configure
+
+ implicit none
+ private
+ save
+
+ !--------------------------------------------------------------------
+ !
+ ! Public parameters
+ !
+ !--------------------------------------------------------------------
+
+ !--------------------------------------------------------------------
+ !
+ ! Public member functions
+ !
+ !--------------------------------------------------------------------
+
+ public :: OcnHorizAdvTracer2ndTend, &
+ OcnHorizAdvTracer2ndInit
+
+ !-------------------------------------------------------------------
+ !
+ ! Private module variables
+ !
+ !--------------------------------------------------------------------
+
+ logical :: &
+ HorizAdv2ndOn !< local flag to determine whether 2nd chosen
+
+ real (kind=RKIND) :: &
+ eddyDiff2 !< base eddy diffusivity for Laplacian
+
+!***********************************************************************
+
+contains
+
+!***********************************************************************
+!
+! routine OcnHorizAdvTracer2ndTend
+!
+!> \brief Computes tendency term for Laplacian horizontal tracer advection
+!> \author Doug Jacobsen
+!> \date 31 August 2011
+!> \version SVN:$Id$
+!> \details
+!> This routine computes the horizontal advection tendency for tracers
+!> based on a Laplacian form for the advection. This tendency takes the
+!> form </font>
<font color="black">abla \cdot (h \kappa_2 </font>
<font color="gray">abla \phi), where \phi is the
+!> tracer, h is the thickness and \kappa_2 is a base horizontal
+!> diffusivity.
+!
+!-----------------------------------------------------------------------
+
+ subroutine OcnHorizAdvTracer2ndTend(tend, tracers, hEdge, u, grid)
+
+ !-----------------------------------------------------------------
+ !
+ ! input variables
+ !
+ !-----------------------------------------------------------------
+
+ real (kind=RKIND), dimension(:,:,:), intent(in) :: &
+ tracers !< Input: array of tracers
+
+ real (kind=RKIND), dimension(:,:), intent(in) :: &
+ hEdge !< Input: layer thickness at cell edges
+
+ real (kind=RKIND), dimension(:,:), intent(in) :: &
+ u !< Input: Velocity
+
+ type (mesh_type), intent(in) :: &
+ grid !< Input: grid information
+
+ !-----------------------------------------------------------------
+ !
+ ! input/output variables
+ !
+ !-----------------------------------------------------------------
+
+ real (kind=RKIND), dimension(:,:,:), intent(inout) :: &
+ tend !< Input/Output: tracer tendency
+
+ !-----------------------------------------------------------------
+ !
+ ! local variables
+ !
+ !-----------------------------------------------------------------
+
+ integer :: &
+ i, k, iCell, iEdge, iTracer, &! loop counters
+ cell1, cell2, &! cell indices
+ nEdges, nCells, nVertLevels, nTracers ! array sizes
+
+ real (kind=RKIND) :: &
+ invAreaCell1, invAreaCell2, &! reciprocal of cell areas
+ flux, tracer_edge ! loop index, flux across edge, and edge tracer value
+
+ real (kind=RKIND), dimension(:), pointer :: &
+ dvEdge, areaCell
+
+ integer, dimension(:), pointer :: maxLevelEdgeTop
+ integer, dimension(:,:), pointer :: cellsOnEdge
+
+ !-----------------------------------------------------------------
+ !
+ ! exit if this advection is not selected
+ !
+ !-----------------------------------------------------------------
+
+ if (.not. HorizAdv2ndOn) return
+
+ !-----------------------------------------------------------------
+ !
+ ! initialize some arrays and sizes
+ !
+ !-----------------------------------------------------------------
+
+ areaCell => grid % areaCell % array
+ cellsOnEdge => grid % cellsOnEdge % array
+ dvEdge => grid % dvEdge % array
+ maxLevelEdgeTop => grid % maxLevelEdgeTop % array
+
+ nTracers = size(tracers,dim=1)
+
+ do iEdge=1,grid % nEdges
+ cell1 = cellsOnEdge(1,iEdge)
+ cell2 = cellsOnEdge(2,iEdge)
+
+ invAreaCell1 = 1.0/areaCell(cell1)
+ invAreaCell2 = 1.0/areaCell(cell2)
+
+ do k=1,maxLevelEdgeTop(iEdge)
+ do iTracer=1,nTracers
+ tracer_edge = 0.5 * (tracers(iTracer,k,cell1) + tracers(iTracer,k,cell2))
+ flux = u(k,iEdge) * dvEdge(iEdge) * hEdge(k,iEdge) * tracer_edge
+ tend(iTracer,k,cell1) = tend(iTracer,k,cell1) - flux*invAreaCell1
+ tend(iTracer,k,cell2) = tend(iTracer,k,cell2) + flux*invAreaCell2
+ end do
+ end do
+ end do
+
+
+ !--------------------------------------------------------------------
+
+ end subroutine OcnHorizAdvTracer2ndTend
+
+!***********************************************************************
+!
+! routine OcnHorizAdvTracer2ndInit
+!
+!> \brief Initializes ocean tracer second order horizontal advection
+!> \author Doug Jacobsen
+!> \date 31 August 2011
+!> \version SVN:$Id$
+!> \details
+!> This routine initializes a variety of quantities related to
+!> second order horizontal tracer advection in the ocean.
+!
+!-----------------------------------------------------------------------
+
+ subroutine OcnHorizAdvTracer2ndInit
+
+ !--------------------------------------------------------------------
+ !
+ ! set some local module variables based on input config choices
+ !
+ !--------------------------------------------------------------------
+
+ HorizAdv2ndOn = .false.
+
+ if (config_tracer_adv_order == 2) then
+ HorizAdv2ndOn = .true.
+ endif
+
+ !--------------------------------------------------------------------
+
+ end subroutine OcnHorizAdvTracer2ndInit
+
+!***********************************************************************
+
+end module OcnHorizAdvTracer2nd
+
+!|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
Added: branches/ocean_projects/performance/mpas/src/core_ocean/module_OcnHorizAdvTracer3rd.F
===================================================================
--- branches/ocean_projects/performance/mpas/src/core_ocean/module_OcnHorizAdvTracer3rd.F (rev 0)
+++ branches/ocean_projects/performance/mpas/src/core_ocean/module_OcnHorizAdvTracer3rd.F 2011-08-31 17:56:45 UTC (rev 969)
@@ -0,0 +1,253 @@
+!|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
+!
+! OcnHorizAdvTracer3rd
+!
+!> \brief Ocean horizontal advection - 3rd order
+!> \author Doug Jacobsen
+!> \date 31 August 2011
+!> \version SVN:$Id:$
+!> \details
+!> This module contains routines for computing horizontal advection
+!> tendencies using a 3rd order formulation.
+!
+!-----------------------------------------------------------------------
+
+module OcnHorizAdvTracer3rd
+
+ use grid_types
+ use configure
+
+ implicit none
+ private
+ save
+
+ !--------------------------------------------------------------------
+ !
+ ! Public parameters
+ !
+ !--------------------------------------------------------------------
+
+ !--------------------------------------------------------------------
+ !
+ ! Public member functions
+ !
+ !--------------------------------------------------------------------
+
+ public :: OcnHorizAdvTracer3rdTend, &
+ OcnHorizAdvTracer3rdInit
+
+ !-------------------------------------------------------------------
+ !
+ ! Private module variables
+ !
+ !--------------------------------------------------------------------
+
+ logical :: &
+ HorizAdv3rdOn !< local flag to determine whether 3rd chosen
+
+ real (kind=RKIND) :: &
+ coef_3rd_order !< coefficient for third order advection
+
+!***********************************************************************
+
+contains
+
+!***********************************************************************
+!
+! routine OcnHorizAdvTracer3rdTend
+!
+!> \brief Computes tendency term for Laplacian horizontal tracer advection
+!> \author Doug Jacobsen
+!> \date 31 August 2011
+!> \version SVN:$Id$
+!> \details
+!> This routine computes the horizontal advection tendency for tracers
+!> based on a Laplacian form for the advection. This tendency takes the
+!> form </font>
<font color="black">abla \cdot (h \kappa_2 </font>
<font color="gray">abla \phi), where \phi is the
+!> tracer, h is the thickness and \kappa_2 is a base horizontal
+!> diffusivity.
+!
+!-----------------------------------------------------------------------
+
+ subroutine OcnHorizAdvTracer3rdTend(tend, tracers, hEdge, u, grid)
+
+ !-----------------------------------------------------------------
+ !
+ ! input variables
+ !
+ !-----------------------------------------------------------------
+
+ real (kind=RKIND), dimension(:,:,:), intent(in) :: &
+ tracers !< Input: array of tracers
+
+ real (kind=RKIND), dimension(:,:), intent(in) :: &
+ hEdge !< Input: layer thickness at cell edges
+
+ real (kind=RKIND), dimension(:,:), intent(in) :: &
+ u !< Input: Velocity
+
+ type (mesh_type), intent(in) :: &
+ grid !< Input: grid information
+
+ !-----------------------------------------------------------------
+ !
+ ! input/output variables
+ !
+ !-----------------------------------------------------------------
+
+ real (kind=RKIND), dimension(:,:,:), intent(inout) :: &
+ tend !< Input/Output: tracer tendency
+
+ !-----------------------------------------------------------------
+ !
+ ! local variables
+ !
+ !-----------------------------------------------------------------
+
+ integer :: &
+ i, k, iCell, iEdge, iTracer, &! loop counters
+ cell1, cell2, &! cell indices
+ nEdges, nCells, nVertLevels, nTracers ! array sizes
+
+ real (kind=RKIND) :: &
+ invAreaCell1, invAreaCell2, &! reciprocal of cell areas
+ flux, tracer_edge, &! flux across edge, and edge tracer value
+ boundaryMask1, boundaryMask2, &
+ d2fdx2_cell1, d2fdx2_cell2
+
+
+ real (kind=RKIND), dimension(:), pointer :: &
+ dcEdge, dvEdge, areaCell
+
+ real (kind=RKIND), dimension(:,:,:), pointer :: deriv_two
+
+ integer, dimension(:), pointer :: maxLevelEdgeTop
+ integer, dimension(:,:), pointer :: cellsOnEdge
+
+
+ !-----------------------------------------------------------------
+ !
+ ! exit if this advection is not selected
+ !
+ !-----------------------------------------------------------------
+
+ if (.not. HorizAdv3rdOn) return
+
+ !-----------------------------------------------------------------
+ !
+ ! initialize some arrays and sizes
+ !
+ !-----------------------------------------------------------------
+
+ areaCell => grid % areaCell % array
+ cellsOnEdge => grid % cellsOnEdge % array
+ dcEdge => grid % dcEdge % array
+ dvEdge => grid % dvEdge % array
+ maxLevelEdgeTop => grid % maxLevelEdgeTop % array
+ deriv_two => grid % deriv_two % array
+
+ nTracers = size(tracers,dim=1)
+
+ do iEdge=1,nEdges
+ cell1 = cellsOnEdge(1,iEdge)
+ cell2 = cellsOnEdge(2,iEdge)
+
+ invAreaCell1 = 1.0 / areaCell(cell1)
+ invAreaCell2 = 1.0 / areaCell(cell2)
+
+ do k=1,maxLevelEdgeTop(iEdge)
+
+ d2fdx2_cell1 = 0.0
+ d2fdx2_cell2 = 0.0
+
+ boundaryMask1 = (.not.grid % boundaryCell % array(k,cell1)) * 1.d0
+ boundaryMask2 = (.not.grid % boundaryCell % array(k,cell2)) * 1.d0
+
+ do iTracer=1,nTracers
+
+ d2fdx2_cell1 = deriv_two(1,1,iEdge) * tracers(iTracer,k,cell1) * boundaryMask1
+ d2fdx2_cell2 = deriv_two(1,2,iEdge) * tracers(iTracer,k,cell2) * boundaryMask2
+
+ !-- all edges of cell 1
+ do i=1, grid % nEdgesOnCell % array (cell1)
+ d2fdx2_cell1 = d2fdx2_cell1 &
+ + deriv_two(i+1,1,iEdge) &
+ * tracers(iTracer,k,grid % CellsOnCell % array (i,cell1)) &
+ * boundaryMask1
+ end do
+
+ !-- all edges of cell 2
+ do i=1, grid % nEdgesOnCell % array (cell2)
+ d2fdx2_cell2 = d2fdx2_cell2 &
+ + deriv_two(i+1,2,iEdge) &
+ * tracers(iTracer,k,grid % CellsOnCell % array (i,cell2)) &
+ * boundaryMask2
+ end do
+
+ !-- if u > 0:
+ if (u(k,iEdge) > 0) then
+ flux = dvEdge(iEdge) * u(k,iEdge) * hEdge(k,iEdge) * ( &
+ 0.5*(tracers(iTracer,k,cell1) + tracers(iTracer,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
+ flux = dvEdge(iEdge) * u(k,iEdge) * hEdge(k,iEdge) * ( &
+ 0.5*(tracers(iTracer,k,cell1) + tracers(iTracer,k,cell2)) &
+ -(dcEdge(iEdge) **2) * (d2fdx2_cell1 + d2fdx2_cell2) / 12. &
+ +(dcEdge(iEdge) **2) * coef_3rd_order*(d2fdx2_cell1 - d2fdx2_cell2) / 12. )
+ end if
+
+ !-- update tendency
+ tend(iTracer,k,cell1) = tend(iTracer,k,cell1) - flux*invAreaCell1
+ tend(iTracer,k,cell2) = tend(iTracer,k,cell2) + flux*invAreaCell2
+ enddo
+ end do
+ end do
+
+ !--------------------------------------------------------------------
+
+ end subroutine OcnHorizAdvTracer3rdTend
+
+!***********************************************************************
+!
+! routine OcnHorizAdvTracer3rdInit
+!
+!> \brief Initializes ocean tracer third order horizontal advection
+!> \author Doug Jacobsen
+!> \date 31 August 2011
+!> \version SVN:$Id$
+!> \details
+!> This routine initializes a variety of quantities related to
+!> third order horizontal tracer advection in the ocean.
+!
+!-----------------------------------------------------------------------
+
+ subroutine OcnHorizAdvTracer3rdInit
+
+ !--------------------------------------------------------------------
+ !
+ ! set some local module variables based on input config choices
+ !
+ !--------------------------------------------------------------------
+
+ HorizAdv3rdOn = .false.
+
+ if (config_tracer_adv_order == 3) then
+ HorizAdv3rdOn = .true.
+ if(config_monotonic) then
+ coef_3rd_order = 0.25
+ else
+ coef_3rd_order = 1.0
+ endif
+ endif
+
+ !--------------------------------------------------------------------
+
+ end subroutine OcnHorizAdvTracer3rdInit
+
+!***********************************************************************
+
+end module OcnHorizAdvTracer3rd
+
+!|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
Added: branches/ocean_projects/performance/mpas/src/core_ocean/module_OcnHorizAdvTracer4th.F
===================================================================
--- branches/ocean_projects/performance/mpas/src/core_ocean/module_OcnHorizAdvTracer4th.F (rev 0)
+++ branches/ocean_projects/performance/mpas/src/core_ocean/module_OcnHorizAdvTracer4th.F 2011-08-31 17:56:45 UTC (rev 969)
@@ -0,0 +1,228 @@
+!|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
+!
+! OcnHorizAdvTracer4th
+!
+!> \brief Ocean horizontal advection - 4th order
+!> \author Doug Jacobsen
+!> \date 31 August 2011
+!> \version SVN:$Id:$
+!> \details
+!> This module contains routines for computing horizontal advection
+!> tendencies using a 4th order formulation.
+!
+!-----------------------------------------------------------------------
+
+module OcnHorizAdvTracer4th
+
+ use grid_types
+ use configure
+
+ implicit none
+ private
+ save
+
+ !--------------------------------------------------------------------
+ !
+ ! Public parameters
+ !
+ !--------------------------------------------------------------------
+
+ !--------------------------------------------------------------------
+ !
+ ! Public member functions
+ !
+ !--------------------------------------------------------------------
+
+ public :: OcnHorizAdvTracer4thTend, &
+ OcnHorizAdvTracer4thInit
+
+ !-------------------------------------------------------------------
+ !
+ ! Private module variables
+ !
+ !--------------------------------------------------------------------
+
+ logical :: &
+ HorizAdv4thOn !< local flag to determine whether 4th chosen
+
+!***********************************************************************
+
+contains
+
+!***********************************************************************
+!
+! routine OcnHorizAdvTracer4thTend
+!
+!> \brief Computes tendency term for Laplacian horizontal tracer advection
+!> \author Doug Jacobsen
+!> \date 31 August 2011
+!> \version SVN:$Id$
+!> \details
+!> This routine computes the horizontal advection tendency for tracers
+!> based on a fourth order form for the advection.
+!
+!-----------------------------------------------------------------------
+
+ subroutine OcnHorizAdvTracer4thTend(tend, tracers, hEdge, u, grid)
+
+ !-----------------------------------------------------------------
+ !
+ ! input variables
+ !
+ !-----------------------------------------------------------------
+
+ real (kind=RKIND), dimension(:,:,:), intent(in) :: &
+ tracers !< Input: array of tracers
+
+ real (kind=RKIND), dimension(:,:), intent(in) :: &
+ hEdge !< Input: layer thickness at cell edges
+
+ real (kind=RKIND), dimension(:,:), intent(in) :: &
+ u !< Input: Velocity
+
+ type (mesh_type), intent(in) :: &
+ grid !< Input: grid information
+
+ !-----------------------------------------------------------------
+ !
+ ! input/output variables
+ !
+ !-----------------------------------------------------------------
+
+ real (kind=RKIND), dimension(:,:,:), intent(inout) :: &
+ tend !< Input/Output: tracer tendency
+
+ !-----------------------------------------------------------------
+ !
+ ! local variables
+ !
+ !-----------------------------------------------------------------
+
+ integer :: &
+ i, k, iCell, iEdge, iTracer, &! loop counters
+ cell1, cell2, &! cell indices
+ nEdges, nCells, nVertLevels, nTracers ! array sizes
+
+ real (kind=RKIND) :: &
+ invAreaCell1, invAreaCell2, &! reciprocal of cell areas
+ flux, tracer_edge, &! flux across edge, and edge tracer value
+ boundaryMask1, boundaryMask2, &
+ d2fdx2_cell1, d2fdx2_cell2
+
+ real (kind=RKIND), dimension(:), pointer :: &
+ dcEdge, dvEdge, areaCell
+
+ real (kind=RKIND), dimension(:,:,:), pointer :: deriv_two
+
+ integer, dimension(:), pointer :: maxLevelEdgeTop
+ integer, dimension(:,:), pointer :: cellsOnEdge
+
+ !-----------------------------------------------------------------
+ !
+ ! exit if this advection is not selected
+ !
+ !-----------------------------------------------------------------
+
+ if (.not. HorizAdv4thOn) return
+
+ !-----------------------------------------------------------------
+ !
+ ! initialize some arrays and sizes
+ !
+ !-----------------------------------------------------------------
+
+ areaCell => grid % areaCell % array
+ cellsOnEdge => grid % cellsOnEdge % array
+ dcEdge => grid % dcEdge % array
+ dvEdge => grid % dvEdge % array
+ maxLevelEdgeTop => grid % maxLevelEdgeTop % array
+ deriv_two => grid % deriv_two % array
+
+ nTracers = size(tracers,dim=1)
+
+ do iEdge=1,nEdges
+ cell1 = cellsOnEdge(1,iEdge)
+ cell2 = cellsOnEdge(2,iEdge)
+
+ invAreaCell1 = 1.0 / areaCell(cell1)
+ invAreaCell2 = 1.0 / areaCell(cell2)
+
+ do k=1,maxLevelEdgeTop(iEdge)
+
+ d2fdx2_cell1 = 0.0
+ d2fdx2_cell2 = 0.0
+ boundaryMask1 = (.not. grid % boundaryCell % array(k,cell1)) * 1.d0
+ boundaryMask2 = (.not. grid % boundaryCell % array(k,cell2)) * 1.d0
+
+ do iTracer=1,nTracers
+ d2fdx2_cell1 = deriv_two(1,1,iEdge) * tracers(iTracer,k,cell1) * boundaryMask1
+ d2fdx2_cell2 = deriv_two(1,2,iEdge) * tracers(iTracer,k,cell2) * boundaryMask2
+
+ !-- all edges of cell 1
+ do i=1, grid % nEdgesOnCell % array (cell1)
+ d2fdx2_cell1 = d2fdx2_cell1 &
+ + deriv_two(i+1,1,iEdge) &
+ * tracers(iTracer,k,grid % CellsOnCell % array (i,cell1)) &
+ * boundaryMask1
+ end do
+
+ !-- all edges of cell 2
+ do i=1, grid % nEdgesOnCell % array (cell2)
+ d2fdx2_cell2 = d2fdx2_cell2 &
+ + deriv_two(i+1,2,iEdge) &
+ * tracers(iTracer,k,grid % CellsOnCell % array (i,cell2)) &
+ * boundaryMask2
+ end do
+
+ flux = dvEdge(iEdge) * u(k,iEdge) * hEdge(k,iEdge) * ( &
+ 0.5*(tracers(iTracer,k,cell1) + tracers(iTracer,k,cell2)) &
+ -(dcEdge(iEdge) **2) * (d2fdx2_cell1 + d2fdx2_cell2) / 12. )
+
+ !-- update tendency
+ tend(iTracer,k,cell1) = tend(iTracer,k,cell1) - flux*invAreaCell1
+ tend(iTracer,k,cell2) = tend(iTracer,k,cell2) + flux*invAreaCell2
+ enddo
+ end do
+ end do
+
+ !--------------------------------------------------------------------
+
+ end subroutine OcnHorizAdvTracer4thTend
+
+!***********************************************************************
+!
+! routine OcnHorizAdvTracer4thInit
+!
+!> \brief Initializes ocean tracer fourth order horizontal advection
+!> \author Doug Jacobsen
+!> \date 31 August 2011
+!> \version SVN:$Id$
+!> \details
+!> This routine initializes a variety of quantities related to
+!> fourth order horizontal tracer advection in the ocean.
+!
+!-----------------------------------------------------------------------
+
+ subroutine OcnHorizAdvTracer4thInit
+
+ !--------------------------------------------------------------------
+ !
+ ! set some local module variables based on input config choices
+ !
+ !--------------------------------------------------------------------
+
+ HorizAdv4thOn = .false.
+
+ if (config_tracer_adv_order == 4) then
+ HorizAdv4thOn = .true.
+ endif
+
+ !--------------------------------------------------------------------
+
+ end subroutine OcnHorizAdvTracer4thInit
+
+!***********************************************************************
+
+end module OcnHorizAdvTracer4th
+
+!|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
Modified: branches/ocean_projects/performance/mpas/src/core_ocean/module_OcnThickness.F
===================================================================
--- branches/ocean_projects/performance/mpas/src/core_ocean/module_OcnThickness.F 2011-08-30 22:18:11 UTC (rev 968)
+++ branches/ocean_projects/performance/mpas/src/core_ocean/module_OcnThickness.F 2011-08-31 17:56:45 UTC (rev 969)
@@ -126,8 +126,8 @@
cell1 = cellsOnEdge(1,iEdge)
cell2 = cellsOnEdge(2,iEdge)
- invAreaCell1 = 1.0/areaCell(cell1)
- invAreaCell2 = 1.0/areaCell(cell2)
+ invAreaCell1 = 1.d0/areaCell(cell1)
+ invAreaCell2 = 1.d0/areaCell(cell2)
do k=1,nVertLevels
flux = u(k,iEdge) * dvEdge(iEdge) * hEdge(k,iEdge)
tend(k,cell1) = tend(k,cell1) - flux * invAreaCell1
@@ -140,8 +140,8 @@
do iEdge=1,nEdges
cell1 = cellsOnEdge(1,iEdge)
cell2 = cellsOnEdge(2,iEdge)
- invAreaCell1 = 1.0/areaCell(cell1)
- invAreaCell2 = 1.0/areaCell(cell2)
+ invAreaCell1 = 1.d0/areaCell(cell1)
+ invAreaCell2 = 1.d0/areaCell(cell2)
do k=1,min(1,maxLevelEdgeTop(iEdge))
flux = u(k,iEdge) * dvEdge(iEdge) * hEdge(k,iEdge)
tend(k,cell1) = tend(k,cell1) - flux * invAreaCell1
@@ -156,7 +156,6 @@
end do
endif ! coordinate type
-
!--------------------------------------------------------------------
end subroutine OcnThicknessTend
Modified: branches/ocean_projects/performance/mpas/src/core_ocean/module_OcnVelocityForcing.F
===================================================================
--- branches/ocean_projects/performance/mpas/src/core_ocean/module_OcnVelocityForcing.F 2011-08-30 22:18:11 UTC (rev 968)
+++ branches/ocean_projects/performance/mpas/src/core_ocean/module_OcnVelocityForcing.F 2011-08-31 17:56:45 UTC (rev 969)
@@ -148,9 +148,14 @@
!-------------------------------------------------------
k = max(maxLevelEdgeTop(iEdge), 1)
+ ! wind stress forcing in top layer only
tend(1,iEdge) = tend(1, iEdge) &
+ u_src(1,iEdge)/rho_ref/hEdge(1,iEdge)
+
+ ! bottom drag is the same as POP:
+ ! -c |u| u where c is unitless and 1.0e-3.
+ ! see POP Reference guide, section 3.4.4.
tend(k,iedge) = tend(k,iEdge) &
- 1.0e-3 * u(k,iEdge) &
* sqrt(2.0 * keEdge(k,iEdge))/hEdge(k,iEdge)
Modified: branches/ocean_projects/performance/mpas/src/core_ocean/module_OcnVmixTracer.F
===================================================================
--- branches/ocean_projects/performance/mpas/src/core_ocean/module_OcnVmixTracer.F 2011-08-30 22:18:11 UTC (rev 968)
+++ branches/ocean_projects/performance/mpas/src/core_ocean/module_OcnVmixTracer.F 2011-08-31 17:56:45 UTC (rev 969)
@@ -4,7 +4,7 @@
!
!> \brief MPAS ocean vertical tracer mixing driver
!> \author Doug Jacobsen
-!> \date 26 July 2011
+!> \date 26 August 2011
!> \version SVN:$Id:$
!> \details
!> This module contains the main driver routine for computing
@@ -58,7 +58,7 @@
!
!> \brief Computes tendency term for vertical tracer mixing
!> \author Doug Jacobsen
-!> \date 26 July 2011
+!> \date 26 August 2011
!> \version SVN:$Id$
!> \details
!> This routine computes the vertical mixing tendency for tracers
@@ -129,7 +129,7 @@
!
!> \brief Initializes ocean tracer vertical mixing quantities
!> \author Doug Jacobsen
-!> \date 26 July 2011
+!> \date 26 August 2011
!> \version SVN:$Id$
!> \details
!> This routine initializes a variety of quantities related to
Modified: branches/ocean_projects/performance/mpas/src/core_ocean/module_OcnVmixTracerConst.F
===================================================================
--- branches/ocean_projects/performance/mpas/src/core_ocean/module_OcnVmixTracerConst.F 2011-08-30 22:18:11 UTC (rev 968)
+++ branches/ocean_projects/performance/mpas/src/core_ocean/module_OcnVmixTracerConst.F 2011-08-31 17:56:45 UTC (rev 969)
@@ -4,7 +4,7 @@
!
!> \brief Ocean vertical mixing - constant parameterization
!> \author Doug Jacobsen
-!> \date 1 August 2011
+!> \date 26 August 2011
!> \version SVN:$Id:$
!> \details
!> This module contains routines for computing vertical mixing
@@ -58,7 +58,7 @@
!
!> \brief Computes tendency term for constant vertical tracer mixing
!> \author Doug Jacobsen
-!> \date 1 August 2011
+!> \date 26 August 2011
!> \version SVN:$Id$
!> \details
!> This routine computes the vertical mixing tendency for tracers
Modified: branches/ocean_projects/performance/mpas/src/core_ocean/module_OcnVmixTracerTanh.F
===================================================================
--- branches/ocean_projects/performance/mpas/src/core_ocean/module_OcnVmixTracerTanh.F 2011-08-30 22:18:11 UTC (rev 968)
+++ branches/ocean_projects/performance/mpas/src/core_ocean/module_OcnVmixTracerTanh.F 2011-08-31 17:56:45 UTC (rev 969)
@@ -4,7 +4,7 @@
!
!> \brief Ocean vertical mixing - hyperbolic tangent parameterization
!> \author Doug Jacobsen
-!> \date 1 August 2011
+!> \date 26 August 2011
!> \version SVN:$Id:$
!> \details
!> This module contains routines for computing vertical mixing
@@ -58,7 +58,7 @@
!
!> \brief Computes tendency term for hyperbolic tangent vertical tracer mixing
!> \author Doug Jacobsen
-!> \date 1 August 2011
+!> \date 26 August 2011
!> \version SVN:$Id$
!> \details
!> This routine computes the vertical mixing tendency for tracers
@@ -138,8 +138,6 @@
maxLevelCell => grid % maxLevelCell % array
zTopZLevel => grid % zTopZLevel % array
- nCellsSolve = grid % nCellsSolve
-
nTracers = size(tracers,dim=1)
nVertLevels = size(tracers, dim=2)
@@ -154,7 +152,7 @@
enddo
fluxVertTop(:,1) = 0.0
- do iCell=1,nCellsSolve
+ do iCell=1,grid % nCellsSolve
do k=2,maxLevelCell(iCell)
do iTracer=1,nTracers
! compute \kappa_v
Modified: branches/ocean_projects/performance/mpas/src/core_ocean/module_mpas_core.F
===================================================================
--- branches/ocean_projects/performance/mpas/src/core_ocean/module_mpas_core.F 2011-08-30 22:18:11 UTC (rev 968)
+++ branches/ocean_projects/performance/mpas/src/core_ocean/module_mpas_core.F 2011-08-31 17:56:45 UTC (rev 969)
@@ -8,6 +8,7 @@
use OcnVmixVel
use OcnVmixTracer
use OcnVelVertAdv
+ use OcnHorizAdvTracer
type (io_output_object) :: restart_obj
integer :: restart_frame
@@ -87,6 +88,7 @@
call OcnVmixVelInit
call OcnVelVertAdvInit
+ call OcnHorizAdvTracerInit
! mrp 100316 In order for this to work, we need to pass domain % dminfo as an
! input arguement into mpas_init. Ask about that later. For now, there will be
Modified: branches/ocean_projects/performance/mpas/src/core_ocean/module_time_integration.F
===================================================================
--- branches/ocean_projects/performance/mpas/src/core_ocean/module_time_integration.F 2011-08-30 22:18:11 UTC (rev 968)
+++ branches/ocean_projects/performance/mpas/src/core_ocean/module_time_integration.F 2011-08-31 17:56:45 UTC (rev 969)
@@ -18,6 +18,7 @@
use OcnThickness
use OcnVelocityForcing
use OcnVelVertAdv
+ use OcnHorizAdvTracer
contains
@@ -443,135 +444,8 @@
! and then change maxLevelEdgeTop to maxLevelEdgeBot in the following section.
! tracer_edge at the boundary will also need to be defined for flux boundaries.
- coef_3rd_order = 0.
- if (config_tracer_adv_order == 3) coef_3rd_order = 1.0
- if (config_tracer_adv_order == 3 .and. config_monotonic) coef_3rd_order = 0.25
+ call OcnHorizAdvTracerTend(tend_tr, tracers, h_edge, u, grid)
- if (config_tracer_adv_order == 2) then
-
- do iEdge=1,nEdges
- cell1 = cellsOnEdge(1,iEdge)
- cell2 = cellsOnEdge(2,iEdge)
-
- invAreaCell1 = 1.0/areaCell(cell1)
- invAreaCell2 = 1.0/areaCell(cell2)
-
- do k=1,maxLevelEdgeTop(iEdge)
- do iTracer=1,num_tracers
- tracer_edge = 0.5 * (tracers(iTracer,k,cell1) + tracers(iTracer,k,cell2))
- flux = u(k,iEdge) * dvEdge(iEdge) * h_edge(k,iEdge) * tracer_edge
- tend_tr(iTracer,k,cell1) = tend_tr(iTracer,k,cell1) - flux*invAreaCell1
- tend_tr(iTracer,k,cell2) = tend_tr(iTracer,k,cell2) + flux*invAreaCell2
- end do
- end do
- end do
-
- else if (config_tracer_adv_order == 3) then
-
- do iEdge=1,nEdges
- cell1 = cellsOnEdge(1,iEdge)
- cell2 = cellsOnEdge(2,iEdge)
-
- invAreaCell1 = 1.0 / areaCell(cell1)
- invAreaCell2 = 1.0 / areaCell(cell2)
-
- do k=1,maxLevelEdgeTop(iEdge)
-
- d2fdx2_cell1 = 0.0
- d2fdx2_cell2 = 0.0
-
- do iTracer=1,num_tracers
-
- !-- 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) * tracers(iTracer,k,cell1)
- d2fdx2_cell2 = deriv_two(1,2,iEdge) * tracers(iTracer,k,cell2)
-
- !-- all edges of cell 1
- do i=1, grid % nEdgesOnCell % array (cell1)
- d2fdx2_cell1 = d2fdx2_cell1 + &
- deriv_two(i+1,1,iEdge) * tracers(iTracer,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) * tracers(iTracer,k,grid % CellsOnCell % array (i,cell2))
- end do
-
- endif
-
- !-- if u > 0:
- if (u(k,iEdge) > 0) then
- flux = dvEdge(iEdge) * u(k,iEdge) * h_edge(k,iEdge) * ( &
- 0.5*(tracers(iTracer,k,cell1) + tracers(iTracer,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
- flux = dvEdge(iEdge) * u(k,iEdge) * h_edge(k,iEdge) * ( &
- 0.5*(tracers(iTracer,k,cell1) + tracers(iTracer,k,cell2)) &
- -(dcEdge(iEdge) **2) * (d2fdx2_cell1 + d2fdx2_cell2) / 12. &
- +(dcEdge(iEdge) **2) * coef_3rd_order*(d2fdx2_cell1 - d2fdx2_cell2) / 12. )
- end if
-
- !-- update tendency
- tend_tr(iTracer,k,cell1) = tend_tr(iTracer,k,cell1) - flux*invAreaCell1
- tend_tr(iTracer,k,cell2) = tend_tr(iTracer,k,cell2) + flux*invAreaCell2
- enddo
- end do
- end do
-
- else if (config_tracer_adv_order == 4) then
-
- do iEdge=1,nEdges
- cell1 = cellsOnEdge(1,iEdge)
- cell2 = cellsOnEdge(2,iEdge)
-
- invAreaCell1 = 1.0 / areaCell(cell1)
- invAreaCell2 = 1.0 / areaCell(cell2)
-
- do k=1,maxLevelEdgeTop(iEdge)
-
- d2fdx2_cell1 = 0.0
- d2fdx2_cell2 = 0.0
-
- do iTracer=1,num_tracers
-
- !-- 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) * tracers(iTracer,k,cell1)
- d2fdx2_cell2 = deriv_two(1,2,iEdge) * tracers(iTracer,k,cell2)
-
- !-- all edges of cell 1
- do i=1, grid % nEdgesOnCell % array (cell1)
- d2fdx2_cell1 = d2fdx2_cell1 + &
- deriv_two(i+1,1,iEdge) * tracers(iTracer,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) * tracers(iTracer,k,grid % CellsOnCell % array (i,cell2))
- end do
-
- endif
-
- flux = dvEdge(iEdge) * u(k,iEdge) * h_edge(k,iEdge) * ( &
- 0.5*(tracers(iTracer,k,cell1) + tracers(iTracer,k,cell2)) &
- -(dcEdge(iEdge) **2) * (d2fdx2_cell1 + d2fdx2_cell2) / 12. )
-
- !-- update tendency
- tend_tr(iTracer,k,cell1) = tend_tr(iTracer,k,cell1) - flux*invAreaCell1
- tend_tr(iTracer,k,cell2) = tend_tr(iTracer,k,cell2) + flux*invAreaCell2
- enddo
- end do
- end do
-
- endif ! if (config_tracer_adv_order == 2 )
-
call timer_stop("scalar horz adv")
call timer_start("scalar vert adv")
</font>
</pre>