<p><b>mhoffman@lanl.gov</b> 2013-02-06 16:04:47 -0700 (Wed, 06 Feb 2013)</p><p>BRANCH COMMIT - land ice<br>
<br>
A simple first-order upwind tendency calculation for tracers. Handy because it works on margin advance.<br>
</p><hr noshade><pre><font color="gray">Modified: branches/land_ice_projects/implement_core/src/core_land_ice/mpas_land_ice_tendency.F
===================================================================
--- branches/land_ice_projects/implement_core/src/core_land_ice/mpas_land_ice_tendency.F 2013-02-06 22:13:05 UTC (rev 2449)
+++ branches/land_ice_projects/implement_core/src/core_land_ice/mpas_land_ice_tendency.F 2013-02-06 23:04:47 UTC (rev 2450)
@@ -260,6 +260,26 @@
tracer_tendency = 0.0_RKIND
select case (config_tracer_advection)
+ case ('FO') !===================================================
+ ! Assign pointers, etc.
+ layerThicknessEdge => state % layerThicknessEdge % array
+ normalVelocity => state % normalVelocity % array
+ tracers => state % tracers % array
+ allocate(uh(mesh % nVertLevels, mesh % nEdges + 1))
+
+ ! FO code needs u*h on edges - this could be calculated elsewhere and saved, potentially (e.g. it's already calculated locally in the FO upwind thickness routine)
+ ! layerThicknessEdge and normalVelocity should match those used in thickness advection!
+ do iEdge = 1, mesh % nEdges
+ do k = 1, mesh % nVertLevels
+ uh(k, iEdge) = normalVelocity(k, iEdge) * layerThicknessEdge(k, iEdge)
+ end do
+ end do
+
+ ! Call the FO!
+ call tracer_advection_tend_fo(tracers, uh, mesh, tracer_tendency, err)
+
+ deallocate(uh)
+
case ('FCT') !===================================================
! Assign pointers, etc.
layerThickness => state % layerThickness % array
@@ -396,6 +416,13 @@
beta = 0.0_RKIND
end where
+ ! Give non-ice cells a temperature of 0
+ do iCell = 1, nCells
+ if ( MASK_IS_NOT_ICE(cellMask(iCell)) ) then
+ state % tracers % array( state%index_temperature, :, iCell) = 0.0_RKIND
+ end if
+ end do
+
! Lower surface is based on floatation for floating ice. For grounded ice it is the bed.
where ( MASK_IS_FLOATING(cellMask) )
lowerSurface = config_sealevel - thickness * (config_ice_density / config_ocean_density)
@@ -1005,5 +1032,117 @@
+
+!***********************************************************************
+!
+! subroutine tracer_advection_tend_fo
+!
+!> \brief tracer_advection_tend_fo
+!> \author Matt Hoffman
+!> \date 06 February 2013
+!> \version SVN:$Id$
+!> \details
+!>
+
+!
+!-----------------------------------------------------------------------
+ subroutine tracer_advection_tend_fo(tracers, uh, mesh, tracer_tendency, err)
+ !-----------------------------------------------------------------
+ !
+ ! input variables
+ !
+ !-----------------------------------------------------------------
+
+ type (mesh_type), intent(in) :: &
+ mesh !< Input: grid information
+
+ real (kind=RKIND), dimension(:,:,:), intent(in) :: &
+ tracers !< Input:
+
+ real (kind=RKIND), dimension(:,:), intent(in) :: &
+ uh !< Input:
+
+ !-----------------------------------------------------------------
+ !
+ ! input/output variables
+ !
+ !-----------------------------------------------------------------
+
+ real (kind=RKIND), dimension(:,:,:), intent(inout) :: &
+ tracer_tendency !< Input: layer tracer tendencies
+
+ !-----------------------------------------------------------------
+ !
+ ! output variables
+ !
+ !-----------------------------------------------------------------
+
+ integer, intent(out) :: err !< Output: error flag
+
+ !-----------------------------------------------------------------
+ !
+ ! local variables
+ !
+ !-----------------------------------------------------------------
+
+ ! pointers to mesh arrays
+ integer, dimension(:), pointer :: nEdgesOnCell
+ integer, dimension(:,:), pointer :: edgesOnCell, cellsOnEdge
+ real (kind=RKIND), dimension(:), pointer :: areaCell, dvEdge
+ ! counters, mesh variables, index variables
+ integer :: nTracers, nCells, nVertLevels
+ integer :: iTracer, iCell, iEdge, thisEdge, k
+ ! stuff for making calculations
+ real (kind=RKIND) :: q, flux
+
+ nCells = mesh % nCells
+ nVertLevels = mesh % nVertLevels
+ nTracers = size(tracers, 1)
+
+ nEdgesOnCell => mesh % nEdgesOnCell % array
+ edgesOnCell => mesh % edgesOnCell % array
+ cellsOnEdge => mesh % cellsOnEdge % array
+ areaCell => mesh % areaCell % array
+ dvEdge => mesh % dvEdge % array
+
+ do iTracer = 1, nTracers
+ do iCell = 1, nCells
+ do iEdge = 1, nEdgesOnCell(iCell)
+ ! What is this edge's index?
+ thisEdge = edgesOnCell(iEdge, iCell)
+ do k = 1, nVertLevels
+ ! Find q-edge using in an upwind sense
+ if (cellsOnEdge(1, thisEdge) == iCell ) then
+ if ( uh(k,thisEdge) < 0.0_RKIND) then
+ ! iCell is the downwind cell so use the other cell for q
+ q = tracers(iTracer, k, cellsOnEdge(2, thisEdge) )
+ flux = -1.0_RKIND * dvEdge(thisEdge) * uh(k, thisEdge) * q
+ else
+ ! iCell is up the upwind cell or the flux is 0
+ q = tracers(iTracer, k, iCell)
+ flux = -1.0_RKIND * dvEdge(thisEdge) * uh(k, thisEdge) * q
+ end if
+ else ! if (cellsOnEdge(2, thisEdge) == iCell ) then
+ if ( uh(k,thisEdge) > 0.0_RKIND) then
+ ! iCell is the downwind cell so use the other cell for q
+ q = tracers(iTracer, k, cellsOnEdge(1, thisEdge) )
+ flux = dvEdge(thisEdge) * uh(k, thisEdge) * q
+ else
+ ! iCell is up the upwind cell or the flux is 0
+ q = tracers(iTracer, k, iCell)
+ flux = dvEdge(thisEdge) * uh(k, thisEdge) * q
+ end if
+ end if
+
+ tracer_tendency(iTracer, k, iCell) = tracer_tendency(iTracer, k, iCell) + flux / areaCell(iCell)
+ end do ! vert levs
+ end do ! edges
+ end do ! cells
+ end do ! tracers
+
+ end subroutine tracer_advection_tend_fo
+
+
+
end module land_ice_tendency
</font>
</pre>