<p><b>mhoffman@lanl.gov</b> 2012-06-26 13:52:13 -0600 (Tue, 26 Jun 2012)</p><p>BRANCH COMMIT -- land ice<br>
<br>
Modifications to advection that 1) don't assume that all layers move in the same direction and 2) don't assume that all layers move in the downslope direction.<br>
Also added calculation of lower ice surface for floating ice (which was previously missing).<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 2012-06-25 23:43:39 UTC (rev 2002)
+++ branches/land_ice_projects/implement_core/src/core_land_ice/mpas_land_ice_tendency.F 2012-06-26 19:52:13 UTC (rev 2003)
@@ -13,6 +13,8 @@
!
!-----------------------------------------------------------------------
+#include "mpas_land_ice_mask.inc"
+
module land_ice_tendency
use mpas_grid_types
@@ -342,35 +344,40 @@
! local variables
!
!-----------------------------------------------------------------
- integer :: iCell, iLevel, nCells, nVertLevels, nVertices
+ integer :: iCell, iLevel, nCells, nVertLevels
real (kind=RKIND), dimension(:), pointer :: thickness, upperSurface, &
- lowerSurface, bedTopography, layerThicknessFractions, thicknessEdge
- real (kind=RKIND), dimension(:,:), pointer :: layerThickness, layerThicknessEdge
- integer, dimension(:), pointer :: vertexMask
- integer, dimension(:,:), pointer :: cellsOnVertex, cellsOnEdge
- integer :: keep_vertex, i, k, iEdge, nEdges, cell1, cell2
+ lowerSurface, bedTopography, layerThicknessFractions
+ integer, dimension(:), pointer :: cellMask
+ real (kind=RKIND), dimension(:,:), pointer :: layerThickness
-
nCells = mesh % nCells
- nEdges = mesh % nEdges
- nVertices = mesh % nVertices
nVertLevels = mesh % nVertLevels
layerThicknessFractions => mesh % layerThicknessFractions % array
- cellsOnVertex => mesh % cellsOnVertex % array
- cellsOnEdge => mesh % cellsOnEdge % array
thickness => state % thickness % array
upperSurface => state % upperSurface % array
lowerSurface => state % lowerSurface % array
bedTopography => state % bedTopography % array
layerThickness => state % layerThickness % array
- vertexMask => state % vertexMask % array
- layerThicknessEdge => state % layerThicknessEdge % array
+ cellMask => state % cellMask % array
+ ! Calculate masks - needs to happen before calculating lower surface so we know where the ice is floating
+ call land_ice_calculate_mask(mesh, state, err)
+
! no ice shelves -- lower surface same as bed topography
! Eventually this should consider floatation!
- lowerSurface(:) = bedTopography(:)
+ where ( MASK_IS_FLOATING(cellMask) )
+ lowerSurface = config_sealevel - thickness * (config_ice_density / config_ocean_density)
+ elsewhere
+ lowerSurface = bedTopography
+ end where
+ ! Make sure lowerSurface calculation is reasonable. This check could be deleted once this has been throroughly tested.
+ do iCell = 1, nCells
+ if (lowerSurface(iCell) < bedTopography(iCell)) then
+ write (0,*) 'lowerSurface less than bedTopography at cell:', iCell
+ endif
+ end do
! as always, upper surface is the lower surface plus the thickness
upperSurface(:) = lowerSurface(:) + thickness(:)
@@ -382,7 +389,69 @@
end do
end do
- ! Calculate h_edge. This is used by both thickness and tracer advection.
+ end subroutine land_ice_diagnostic_solve
+
+
+!***********************************************************************
+!
+! subroutine land_ice_diagnostic_solve_using_velocity
+!
+!> \brief Computes diagnostic variables that require knowing velocity
+!> \author Matt Hoffman
+!> \date 19 April 2012
+!> \version SVN:$Id$
+!> \details
+!> This routine computes the diagnostic variables that require knowing velocity for land ice
+!
+!-----------------------------------------------------------------------
+ subroutine land_ice_diagnostic_solve_using_velocity(mesh, state, err)!{{{
+
+ !-----------------------------------------------------------------
+ !
+ ! input variables
+ !
+ !-----------------------------------------------------------------
+
+ type (mesh_type), intent(in) :: &
+ mesh !< Input: grid information
+
+ !-----------------------------------------------------------------
+ !
+ ! input/output variables
+ !
+ !-----------------------------------------------------------------
+
+ type (state_type), intent(inout) :: &
+ state !< Input/Output: state for which to update diagnostic variables
+
+ !-----------------------------------------------------------------
+ !
+ ! output variables
+ !
+ !-----------------------------------------------------------------
+
+ integer, intent(out) :: err !< Output: error flag
+
+ !-----------------------------------------------------------------
+ !
+ ! local variables
+ !
+ !-----------------------------------------------------------------
+ real (kind=RKIND), dimension(:,:), pointer :: layerThickness, layerThicknessEdge, normalVelocity
+ integer, dimension(:,:), pointer :: cellsOnEdge
+ integer :: nVertLevels, iEdge, nEdges, cell1, cell2, k
+ real (kind=RKIND) :: VelSign
+
+ nEdges = mesh % nEdges
+ nVertLevels = mesh % nVertLevels
+ cellsOnEdge => mesh % cellsOnEdge % array
+
+ layerThickness => state % layerThickness % array
+ normalVelocity => state % normalVelocity % array
+ layerThicknessEdge => state % layerThicknessEdge % array
+
+
+ ! Calculate h_edge. This is used by both thickness and tracer advection on the following Forward Euler time step.
! Note: FO-Upwind thickness advection does not explicitly use h_edge but a FO h_edge is implied.
! Note: SIA velocity solver uses its own local calculation of h_edge that is always 2nd order.
! Note: ocn_diagnostic_solve in mpas_ocn_tendency.F has 2, 3, & 4th order calculations for h_edge that can be used.
@@ -396,7 +465,8 @@
cell2 = cellsOnEdge(2,iEdge)
do k=1, nVertLevels
! Calculate h on edges using first order
- layerThicknessEdge(k,iEdge) = max(layerThickness(k, cell1), layerThickness(k, cell2))
+ VelSign = sign(1.0, normalVelocity(k, iEdge))
+ layerThicknessEdge(k,iEdge) = max(VelSign * layerThickness(k, cell1), VelSign * (-1.0) * layerThickness(k, cell2)) ! + velocity goes from index 1 to 2 in the cellsOnEdge array. ! Doug does the calculation as: h_edge = max(VelSign, 0.0) * h1 - min(VelSign, 0.0) * h2
end do
! thickness_edge is not currently in registry and not currenly needed. If it is, uncomment the next line
!h_edge = max(thickness(cell1), thickness(cell2))
@@ -406,139 +476,134 @@
print *,'layerThicknessEdge not calculated!'
endif
- ! Calculate masks
- call land_ice_calculate_mask(mesh, state, err)
+ ! Note: no halo update needed on layerThicknessEdge because normalVelocity has already been halo updated.
+ end subroutine land_ice_diagnostic_solve_using_velocity
- end subroutine land_ice_diagnostic_solve
-
-
-
-
! private subroutines================================================
-!***********************************************************************
-!
-! subroutine land_ice_tend_h_fo_upwind
-!
-!> \brief Computes tendency term from horizontal advection of thickness
-!> \author Matt Hoffman
-!> \date 16 April 2012
-!> \version SVN:$Id$
-!> \details
-!> This routine computes the horizontal advection tendency for
-!> thickness based on current state and user choices of forcings. Based on
-!> ocn_thick_hadv_tend in the ocean core.
-!
-!-----------------------------------------------------------------------
+!!***********************************************************************
+!!
+!! subroutine land_ice_tend_h_fo_upwind
+!!
+!!> \brief Computes tendency term from horizontal advection of thickness
+!!> \author Matt Hoffman
+!!> \date 16 April 2012
+!!> \version SVN:$Id$
+!!> \details
+!!> This routine computes the horizontal advection tendency for
+!!> thickness based on current state and user choices of forcings. Based on
+!!> ocn_thick_hadv_tend in the ocean core.
+!!
+!!-----------------------------------------------------------------------
- subroutine land_ice_tend_h_fo_upwind(grid, normalVelocity, layerThickness, thickness, thickness_tend, dt, err)!{{{
+! subroutine land_ice_tend_h_fo_upwind(grid, normalVelocity, layerThickness, thickness, thickness_tend, dt, err)!{{{
- !-----------------------------------------------------------------
- !
- ! input variables
- !
- !-----------------------------------------------------------------
+! !-----------------------------------------------------------------
+! !
+! ! input variables
+! !
+! !-----------------------------------------------------------------
- real (kind=RKIND), dimension(:,:), intent(in) :: &
- normalVelocity !< Input: velocity
+! real (kind=RKIND), dimension(:,:), intent(in) :: &
+! normalVelocity !< Input: velocity
- real (kind=RKIND), dimension(:,:), intent(in) :: &
- layerThickness !< Input: thickness of each layer
+! real (kind=RKIND), dimension(:,:), intent(in) :: &
+! layerThickness !< Input: thickness of each layer
- real (kind=RKIND), dimension(:), intent(in) :: &
- thickness !< Input: thickness
+! real (kind=RKIND), dimension(:), intent(in) :: &
+! thickness !< Input: thickness
- type (mesh_type), intent(in) :: &
- grid !< Input: grid information
+! type (mesh_type), intent(in) :: &
+! grid !< Input: grid information
- real (kind=RKIND) :: dt
+! real (kind=RKIND) :: dt
- !-----------------------------------------------------------------
- !
- ! input/output variables
- !
- !-----------------------------------------------------------------
+! !-----------------------------------------------------------------
+! !
+! ! input/output variables
+! !
+! !-----------------------------------------------------------------
- real (kind=RKIND), dimension(:), intent(inout) :: &
- thickness_tend !< Input/Output: thickness tendency
+! real (kind=RKIND), dimension(:), intent(inout) :: &
+! thickness_tend !< Input/Output: thickness tendency
- !-----------------------------------------------------------------
- !
- ! output variables
- !
- !-----------------------------------------------------------------
+! !-----------------------------------------------------------------
+! !
+! ! output variables
+! !
+! !-----------------------------------------------------------------
- integer, intent(out) :: err !< Output: error flag
+! integer, intent(out) :: err !< Output: error flag
- !-----------------------------------------------------------------
- !
- ! local variables
- !
- !-----------------------------------------------------------------
+! !-----------------------------------------------------------------
+! !
+! ! local variables
+! !
+! !-----------------------------------------------------------------
- integer :: iEdge, nEdges, cell1, cell2, nVertLevels, k, CellUpwind, CellDownwind
+! integer :: iEdge, nEdges, cell1, cell2, nVertLevels, k, CellUpwind, CellDownwind
- integer, dimension(:,:), pointer :: cellsOnEdge
+! integer, dimension(:,:), pointer :: cellsOnEdge
- real (kind=RKIND) :: flux, VelSign, h_edge, ubar, maxAllowableDt
- real (kind=RKIND), dimension(:), pointer :: dvEdge, areaCell, dcEdge
+! real (kind=RKIND) :: flux, VelSign, h_edge, ubar, maxAllowableDt
+! real (kind=RKIND), dimension(:), pointer :: dvEdge, areaCell, dcEdge
- err = 0
- maxAllowableDt = 1.0e36
+! err = 0
+! maxAllowableDt = 1.0e36
- nEdges = grid % nEdges
- nVertLevels = grid % nVertLevels
- cellsOnEdge => grid % cellsOnEdge % array
- dvEdge => grid % dvEdge % array
- dcEdge => grid % dcEdge % array
- areaCell => grid % areaCell % array
+! nEdges = grid % nEdges
+! nVertLevels = grid % nVertLevels
+! cellsOnEdge => grid % cellsOnEdge % array
+! dvEdge => grid % dvEdge % array
+! dcEdge => grid % dcEdge % array
+! areaCell => grid % areaCell % array
- !print *,'Max velocity magn.:', maxval(abs(normalVelocity))
+! !print *,'Max velocity magn.:', maxval(abs(normalVelocity))
- ! Zero the tendency before accumulating
- thickness_tend = 0.0
+! ! Zero the tendency before accumulating
+! thickness_tend = 0.0
- do iEdge=1,nEdges
- cell1 = cellsOnEdge(1,iEdge)
- cell2 = cellsOnEdge(2,iEdge)
+! do iEdge=1,nEdges
+! cell1 = cellsOnEdge(1,iEdge)
+! cell2 = cellsOnEdge(2,iEdge)
- VelSign = sign(1.0, normalVelocity(1, iEdge))
- if (VelSign .gt. 0.0) then
- CellUpwind = cell1
- CellDownwind = cell2
- else
- CellUpwind = cell2
- CellDownwind = cell1
- endif
- if (thickness(cellUpwind) .gt. 0.0) then ! Don't calculate for non-ice cells - would result in divide by 0
- h_edge = thickness(CellUpwind)
- flux = 0.0
- do k=1, nVertLevels
- ! Calculate thickness averaged velocity - this make be calculated externally and passed in
- flux = flux + layerThickness(k, cellUpwind) * abs(normalVelocity(k, iEdge))
- enddo
- ubar = flux / thickness(cellUpwind)
- if ( (abs(ubar) * dt/SecondsInYear) .gt. (0.5 * dcEdge(iEdge))) then
- !maxAllowableDt = min(maxAllowableDt, (0.5 * dcEdge)/ubar )
- write(0,*) 'CFL violation at edge', iEdge
- err = err + 1
- endif
- thickness_tend(cellUpwind) = thickness_tend(cellUpwind) - flux * dvEdge(iEdge) / areaCell(cellUpwind)
- thickness_tend(cellDownwind) = thickness_tend(cellDownwind) + flux * dvEdge(iEdge) / areaCell(cellDownwind)
- endif
- end do
+! VelSign = sign(1.0, normalVelocity(1, iEdge))
+! if (VelSign .gt. 0.0) then
+! CellUpwind = cell1
+! CellDownwind = cell2
+! else
+! CellUpwind = cell2
+! CellDownwind = cell1
+! endif
+! if (thickness(cellUpwind) .gt. 0.0) then ! Don't calculate for non-ice cells - would result in divide by 0
+! h_edge = thickness(CellUpwind)
+! flux = 0.0
+! do k=1, nVertLevels
+! ! Calculate thickness averaged velocity - this make be calculated externally and passed in
+! flux = flux + layerThickness(k, cellUpwind) * abs(normalVelocity(k, iEdge))
+! enddo
+! ubar = flux / thickness(cellUpwind)
+! if ( (abs(ubar) * dt/SecondsInYear) .gt. (0.5 * dcEdge(iEdge))) then
+! !maxAllowableDt = min(maxAllowableDt, (0.5 * dcEdge)/ubar )
+! write(0,*) 'CFL violation at edge', iEdge
+! err = err + 1
+! endif
+! thickness_tend(cellUpwind) = thickness_tend(cellUpwind) - flux * dvEdge(iEdge) / areaCell(cellUpwind)
+! thickness_tend(cellDownwind) = thickness_tend(cellDownwind) + flux * dvEdge(iEdge) / areaCell(cellDownwind)
+! endif
+! end do
- if (err .gt. 0) then
- write(6,*) 'CFL violation at ', err, ' edges! Maximum time step should be ', maxAllowableDt
- err = 1
- endif
- !--------------------------------------------------------------------
+! if (err .gt. 0) then
+! write(6,*) 'CFL violation at ', err, ' edges! Maximum time step should be ', maxAllowableDt
+! err = 1
+! endif
+! !--------------------------------------------------------------------
- end subroutine land_ice_tend_h_fo_upwind !}}}
+! end subroutine land_ice_tend_h_fo_upwind !}}}
@@ -606,7 +671,7 @@
integer, dimension(:,:), pointer :: cellsOnEdge
- real (kind=RKIND) :: flux, VelSign, maxAllowableDt
+ real (kind=RKIND) :: flux, VelSign, maxAllowableDt, VelSignSum
real (kind=RKIND), dimension(:), pointer :: dvEdge, areaCell, dcEdge
err = 0
@@ -623,16 +688,20 @@
do iEdge=1,nEdges
cell1 = cellsOnEdge(1,iEdge)
cell2 = cellsOnEdge(2,iEdge)
- VelSign = sign(1.0, normalVelocity(1, iEdge))
- if (VelSign .gt. 0.0) then
- CellUpwind = cell1
- CellDownwind = cell2
- else
- CellUpwind = cell2
- CellDownwind = cell1
- endif
- if (layerThickness(1,cellUpwind) .gt. 0.0) then ! Don't calculate for non-ice cells - would result in divide by 0 (just need to check one layer)
+ if ( (layerThickness(1,cell1) .gt. 0.0) .or. (layerThickness(1,cell2) .gt. 0.0) ) then ! Don't calculate for non-ice cells - would result in divide by 0 (just need to check one layer) \todo: this should use the edgeMask and use the dynamic thickness limit
+
+ VelSignSum = 0.0
do k=1, nVertLevels
+ ! Don't assume all layers are flowing the same direction, so determine the upwind direction separately for each layer
+ VelSign = sign(1.0, normalVelocity(k, iEdge))
+ VelSignSum = VelSignSum + VelSign
+ if (VelSign .gt. 0.0) then
+ CellUpwind = cell1
+ CellDownwind = cell2
+ else
+ CellUpwind = cell2
+ CellDownwind = cell1
+ endif
maxAllowableDt = (0.5 * dcEdge(iEdge)) / (abs(normalVelocity(k, iEdge)) + 1.0e-36) ! in years
if ( (dt/SecondsInYear) .gt. maxAllowableDt ) then
!write(0,*) 'CFL violation at level, edge', k, iEdge
@@ -644,6 +713,10 @@
tend(k,cell1) = tend(k,cell1) - flux / areaCell(cell1)
tend(k,cell2) = tend(k,cell2) + flux / areaCell(cell2)
end do
+ ! Sanity check on ice motion:
+ if (abs(VelSignSum) .ne. nVertLevels) then
+ write (6,*) 'Unusual occurrence: Not all layers are moving the same direction on edge ', iEdge
+ endif
end if
end do
Modified: branches/land_ice_projects/implement_core/src/core_land_ice/mpas_land_ice_time_integration_forwardeuler.F
===================================================================
--- branches/land_ice_projects/implement_core/src/core_land_ice/mpas_land_ice_time_integration_forwardeuler.F 2012-06-25 23:43:39 UTC (rev 2002)
+++ branches/land_ice_projects/implement_core/src/core_land_ice/mpas_land_ice_time_integration_forwardeuler.F 2012-06-26 19:52:13 UTC (rev 2003)
@@ -20,6 +20,7 @@
use mpas_constants
use mpas_dmpar
use mpas_timer
+ use mpas_vector_reconstruction
use land_ice_vel, only: land_ice_vel_solve
use land_ice_tendency
@@ -197,7 +198,7 @@
thicknessNew = 0.0
end where
if (sum(mask) .gt. 0) then
- print *,' Cells with negative thickness (set to 0):',sum(mask)
+ print *,' Cells with negative thickness (set to 0):',sum(mask)
endif
mask = 0
@@ -271,16 +272,34 @@
! Compute the diagnostic velocity for this time step using the updated (new) state
! Pass in the old velocity to be used as an initial guess for velocity solvers that need one.
normalVelocityNew = normalVelocityOld
- call land_ice_vel_solve(mesh, stateNew, err)
+ call land_ice_vel_solve(mesh, stateNew, err) ! ****** Calculate Velocity ******
! update halos on velocity
call mpas_dmpar_exch_halo_field2d_real(dminfo, normalVelocityNew, &
nVertLevels, mesh % nEdges, &
block % parinfo % edgesToSend, block % parinfo % edgesToRecv)
+ ! Calculate reconstructed velocities for this time step - do this after velocity halo update in case velocities on the 1-halo edge are wrong (depends on velocity solver)
+ call mpas_reconstruct(mesh, stateNew % normalVelocity % array,&
+ stateNew % uReconstructX % array, &
+ stateNew % uReconstructY % array, &
+ stateNew % uReconstructZ % array, &
+ stateNew % uReconstructZonal % array, &
+ stateNew % uReconstructMeridional % array &
+ )
+
block => block % next
end do
call mpas_timer_stop("velocity solve")
+
+ call mpas_timer_start("calc. diagnostic vars except vel")
+ block => domain % blocklist
+ do while (associated(block))
+ call land_ice_diagnostic_solve(mesh, stateNew, err) ! Some diagnostic variables require velocity to compute
+ block => block % next
+ end do
+ call mpas_timer_stop("calc. diagnostic vars except vel")
+
! === Cleanup & Misc. =============================
Modified: branches/land_ice_projects/implement_core/src/core_land_ice/mpas_land_ice_vel.F
===================================================================
--- branches/land_ice_projects/implement_core/src/core_land_ice/mpas_land_ice_vel.F 2012-06-25 23:43:39 UTC (rev 2002)
+++ branches/land_ice_projects/implement_core/src/core_land_ice/mpas_land_ice_vel.F 2012-06-26 19:52:13 UTC (rev 2003)
@@ -254,7 +254,6 @@
subroutine land_ice_vel_solve(mesh, state, err)
use land_ice_SIA
- use mpas_vector_reconstruction
!-----------------------------------------------------------------
!
@@ -305,14 +304,6 @@
return
end select
- ! Calculate reconstructed velocities for this time step
- call mpas_reconstruct(mesh, state % normalVelocity % array,&
- state % uReconstructX % array, &
- state % uReconstructY % array, &
- state % uReconstructZ % array, &
- state % uReconstructZonal % array, &
- state % uReconstructMeridional % array &
- )
!--------------------------------------------------------------------
</font>
</pre>