<p><b>mhoffman@lanl.gov</b> 2012-05-21 15:14:18 -0600 (Mon, 21 May 2012)</p><p>BRANCH COMMIT -- land ice<br>
<br>
Updated shallow ice velocity solver to use correct vertical indexing convention.<br>
Also cleaned up various implementations of thickness on edges:<br>
SIA velocity is hardcoded to use 2nd order h_edge.<br>
FO-Upwind thickness advection is hardcoded to use 1st order h_edge.<br>
Diagnostic h_edge (needed for tracer advection) is determined by choice of thickness advection scheme.<br>
Also fixed a bug in the Makefile.<br>
</p><hr noshade><pre><font color="gray">Modified: branches/land_ice_projects/implement_core/src/core_land_ice/Makefile
===================================================================
--- branches/land_ice_projects/implement_core/src/core_land_ice/Makefile 2012-05-21 20:07:29 UTC (rev 1924)
+++ branches/land_ice_projects/implement_core/src/core_land_ice/Makefile 2012-05-21 21:14:18 UTC (rev 1925)
@@ -34,13 +34,13 @@
mpas_land_ice_time_integration.o: mpas_land_ice_time_integration_forwardeuler.o
-mpas_land_ice_time_integration_forwardeuler.o: mpas_land_ice_vel.o mpas_land_ice_tendency.o mpas_ocn_tracer_advection.o
+mpas_land_ice_time_integration_forwardeuler.o: mpas_land_ice_vel.o mpas_land_ice_tendency.o
mpas_land_ice_mask.o:
mpas_land_ice_vel.o: mpas_land_ice_lifev.o mpas_land_ice_sia.o
-mpas_land_ice_tendency.o: mpas_land_ice_mask.o
+mpas_land_ice_tendency.o: mpas_land_ice_mask.o mpas_ocn_tracer_advection.o
mpas_land_ice_lifev.o:
Modified: branches/land_ice_projects/implement_core/src/core_land_ice/Registry
===================================================================
--- branches/land_ice_projects/implement_core/src/core_land_ice/Registry 2012-05-21 20:07:29 UTC (rev 1924)
+++ branches/land_ice_projects/implement_core/src/core_land_ice/Registry 2012-05-21 21:14:18 UTC (rev 1925)
@@ -18,7 +18,7 @@
% valid dycores are SIA, L1L2, FO, Stokes. All but SIA require compiling with LifeV libraries.
namelist character land_ice_model config_time_integration ForwardEuler
% valid time integration is ForwardEuler
-namelist character land_ice_model config_advection FO-Upwind
+namelist character land_ice_model config_thickness_advection FO-Upwind
% valid advection is FO-Upwind, None (FCT coming soon)
namelist character land_ice_model config_tracer_advection FCT
% valid options for tracer_advection are FCT (flux-corrected transport), None
@@ -175,8 +175,8 @@
% Mesh variables specific to land ice core
% Sigma coordinate: read from input, saved in restart and written to output
var persistent real layerThicknessFractions ( nVertLevels ) 0 iro layerThicknessFractions mesh - -
-var persistent real layerCenterSigma ( nVertLevels ) 0 - layerCenterSigma mesh - -
-var persistent real layerInterfaceSigma ( nVertLevelsPlus1 ) 0 - layerInterfaceSigma mesh - -
+var persistent real layerCenterSigma ( nVertLevels ) 0 o layerCenterSigma mesh - -
+var persistent real layerInterfaceSigma ( nVertLevelsPlus1 ) 0 o layerInterfaceSigma mesh - -
% Boundary Conditions (defined as mesh variables because they are read in once and held constant for the simulation - at least for now)
%
@@ -199,7 +199,6 @@
% Tendency variables
var persistent real tend_layerThickness ( nVertLevels nCells Time ) 1 o layerThickness tend - -
-var persistent real tend_thickness ( nCells Time ) 1 o thickness tend - -
var persistent real tend_temperature ( nVertLevels nCells Time ) 1 o temperature tend tracers dynamics
%var persistent real tend_sup_thickness ( nVertLevels nCells Time ) 1 o sup_thickness tend sup_thickness bob
Modified: branches/land_ice_projects/implement_core/src/core_land_ice/mpas_land_ice_global_diagnostics.F
===================================================================
--- branches/land_ice_projects/implement_core/src/core_land_ice/mpas_land_ice_global_diagnostics.F 2012-05-21 20:07:29 UTC (rev 1924)
+++ branches/land_ice_projects/implement_core/src/core_land_ice/mpas_land_ice_global_diagnostics.F 2012-05-21 21:14:18 UTC (rev 1925)
@@ -50,7 +50,7 @@
real (kind=RKIND), dimension(:), pointer :: areaCell, dcEdge, dvEdge, areaTriangle
!!! real (kind=RKIND), dimension(:), pointer :: h_s, fCell, fEdge
!!! real (kind=RKIND), dimension(:,:), pointer :: h, u, v, layerThicknessEdge, &
- real (kind=RKIND), dimension(:,:), pointer :: layerThickness, normalVelocity, layerThicknessEdge, &
+ real (kind=RKIND), dimension(:,:), pointer :: layerThickness, normalVelocity, &
!!! pv_edge, pv_vertex, pv_cell, &
layerThicknessVertex, weightsOnEdge
@@ -110,7 +110,7 @@
normalVelocity => state % normalVelocity % array
!!! v => state % v % array
tracers => state % tracers % array
- layerThicknessEdge => state % layerThicknessEdge % array
+! layerThicknessEdge => state % layerThicknessEdge % array
layerThicknessVertex => state % layerThicknessVertex % array
!!! pv_edge => state % pv_edge % array
!!! pv_vertex => state % pv_vertex % array
Modified: branches/land_ice_projects/implement_core/src/core_land_ice/mpas_land_ice_sia.F
===================================================================
--- branches/land_ice_projects/implement_core/src/core_land_ice/mpas_land_ice_sia.F 2012-05-21 20:07:29 UTC (rev 1924)
+++ branches/land_ice_projects/implement_core/src/core_land_ice/mpas_land_ice_sia.F 2012-05-21 21:14:18 UTC (rev 1925)
@@ -215,7 +215,19 @@
end subroutine land_ice_sia_block_init
!***********************************************************************
-
+!
+! subroutine land_ice_sia_solve
+!
+!> \brief Computes velocity using Shallow Ice Appoximation
+!> \author Matt Hoffman
+!> \date 21 May 2012
+!> \version SVN:$Id$
+!> \details
+!> This routine computes the normal velocity on edges for each layer
+!> using the Shallow Ice Approximation. It calculates ice thickness on
+!> on an edge using the average of the two neighoring cells (2nd order).
+!
+!-----------------------------------------------------------------------
subroutine land_ice_sia_solve(mesh, state, err)
use mpas_constants, only: gravity
@@ -251,63 +263,57 @@
!
!-----------------------------------------------------------------
- real (kind=RKIND), dimension(:), pointer :: thickness, layerThicknessFractions
- real (kind=RKIND), dimension(:), pointer :: dcEdge
- real (kind=RKIND), dimension(:,:), pointer :: layerThickness, normalVelocity
+ real (kind=RKIND), dimension(:), pointer :: thickness, layerCenterSigma, dcEdge
+ real (kind=RKIND), dimension(:,:), pointer :: normalVelocity
integer, dimension(:,:), pointer :: cellsOnEdge
integer, dimension(:), pointer :: edgeMask
- real (kind=RKIND), dimension(mesh%nVertLevels) :: layerFractionAboveBed
- integer :: nVertLevels, nCells, nEdges, iCell, iLevel, iEdge
- real (kind=RKIND) :: ratefactor, basalVelocity, slopeOnEdge, thicknessOnEdge
- real rhoi !!!! THIS SHOULD BE A PHYSICAL PARAMETER ELSEWHERE
+ integer :: nVertLevels, nEdges, iLevel, iEdge, cell1, cell2
+ real (kind=RKIND) :: ratefactor, basalVelocity, slopeOnEdge, &
+ layerCenterHeightOnEdge, thicknessEdge
+ real rhoi
err = 0
! Set needed variables and pointers
- nCells = mesh % nCells
nEdges = mesh % nEdges
nVertLevels = mesh % nVertLevels
dcEdge => mesh % dcEdge % array
cellsOnEdge => mesh % cellsOnEdge % array
- layerThicknessFractions => mesh % layerThicknessFractions % array
+ layerCenterSigma => mesh % layerCenterSigma % array
normalVelocity => state % normalVelocity % array
thickness => state % thickness % array
- layerThickness => state % layerThickness % array
edgeMask => state % edgeMask % array
+
rhoi = config_ice_density
basalVelocity = 0.0 ! Assume no sliding
+ ! Calculate ratefactor (A) at edge - This should be calculated external to this subroutine and as a function of temperature
ratefactor = 1.0e-16
- ! Determine the height of each layer above the bed
- layerFractionAboveBed(1) = 0.5 * layerThicknessFractions(1)
- do iLevel = 2, nVertLevels
- layerFractionAboveBed(iLevel) = layerFractionAboveBed(iLevel - 1) + &
- 0.5 * layerThicknessFractions(iLevel - 1) + 0.5 * layerThicknessFractions(iLevel)
- end do
-
! Loop over edges
do iEdge = 1, nEdges
! Only calculate velocity for edges that are part of the ice sheet.
! Also, the velocity calculation should be valid for non-ice edges (i.e. returns 0).
if ( MASK_IS_ICE(edgeMask(iEdge)) ) then
- ! Calculate slope, thickness at edge
- ! These could/should be calculated externally to this subroutine
- slopeOnEdge = (thickness(cellsOnEdge(1,iEdge)) - thickness(cellsOnEdge(2,iEdge)) ) &
- / dcEdge(iEdge)
- thicknessOnEdge = 0.5 * (thickness(cellsOnEdge(1,iEdge)) + thickness(cellsOnEdge(2,iEdge)))
+ cell1 = cellsOnEdge(1,iEdge)
+ cell2 = cellsOnEdge(2,iEdge)
+ ! Calculate slope at edge
+ ! This could/should be calculated externally to this subroutine
+ slopeOnEdge = (thickness(cell1) - thickness(cell2) ) / dcEdge(iEdge)
+ ! Calculate thickness on edge - 2nd order
+ thicknessEdge = (thickness(cell1) + thickness(cell2) ) * 0.5
! Loop over layers
do iLevel = 1, nVertLevels
- ! Calculate ratefactor (A) at edge - This should be calculated here as a function of temperature
- ! ratefactor = f(T)
+ ! Determine the height of each layer above the bed
+ layerCenterHeightOnEdge = thicknessEdge * (1.0 - layerCenterSigma(iLevel) )
! Calculate SIA velocity
normalVelocity(iLevel,iEdge) = basalVelocity + &
0.5 * ratefactor * (rhoi * gravity)**3 * slopeOnEdge**3 * &
- (thicknessOnEdge**4 - (thicknessOnEdge - thicknessOnEdge * layerFractionAboveBed(iLevel))**4)
+ (thicknessEdge**4 - (thicknessEdge - layerCenterHeightOnEdge)**4)
end do ! Levels
endif
end do ! edges
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-05-21 20:07:29 UTC (rev 1924)
+++ branches/land_ice_projects/implement_core/src/core_land_ice/mpas_land_ice_tendency.F 2012-05-21 21:14:18 UTC (rev 1925)
@@ -113,7 +113,7 @@
! 0 tendency
layerThickness_tend = 0.0
- select case (config_advection)
+ select case (config_thickness_advection)
case ('FO-Upwind') !===================================================
!print *,'Using FO Upwind for thickness advection'
@@ -122,7 +122,7 @@
! state % layerThickness % array, state % thickness % array, thickness_tend, dt, err)
call land_ice_tend_h_layers_fo_upwind(mesh, state % normalVelocity % array, &
- state % layerThicknessEdge % array, layerThickness_tend, dt, err)
+ state % layerThickness % array, layerThickness_tend, dt, err)
!case ('FCT') !===================================================
! MJH TEMP TRYING IT WITH THICKNESS
@@ -142,7 +142,7 @@
case ('None') !===================================================
! Do nothing
case default !===================================================
- write(0,*) trim(config_advection), ' is not a valid thickness advection option.'
+ write(0,*) trim(config_thickness_advection), ' is not a valid thickness advection option.'
call mpas_dmpar_abort(dminfo)
end select !===================================================
@@ -246,6 +246,7 @@
end where
! FCT 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)
@@ -332,7 +333,7 @@
!-----------------------------------------------------------------
integer :: iCell, iLevel, nCells, nVertLevels, nVertices
real (kind=RKIND), dimension(:), pointer :: thickness, upperSurface, &
- lowerSurface, bedTopography, layerThicknessFractions
+ lowerSurface, bedTopography, layerThicknessFractions, thicknessEdge
real (kind=RKIND), dimension(:,:), pointer :: layerThickness, layerThicknessEdge
integer, dimension(:), pointer :: vertexMask
integer, dimension(:,:), pointer :: cellsOnVertex, cellsOnEdge
@@ -370,29 +371,31 @@
end do
end do
-
- ! \todo Calculate h_edge that can be used by SIA solve and FO advection scheme.
- ! ocn_diagnostic_solve in mpas_ocn_tendency.F has 2, 3, & 4th order calculations for h_edge that can be used.
- ! Right now I am using 2nd order h_edge in SIA solve and 1st order h_edge in FO advection.
+ ! Calculate h_edge. This is used by both thickness and tracer advection.
+ ! 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.
! NOTE: This calculates FO upwind h edge
! Both thickness and layerThickness should be updated by this time.
- do iEdge=1,nEdges
- cell1 = cellsOnEdge(1,iEdge)
- 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))
+ if (config_thickness_advection .eq. 'FO-Upwind') then
+ ! If using FO-Upwind then h_edge must be FO.
+ do iEdge=1,nEdges
+ cell1 = cellsOnEdge(1,iEdge)
+ 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))
+ 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))
+ !!!h_edge = (thickness(k) + thickness(k) ) / 2.0 ! 2nd order
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))
- !!!h_edge = (thickness(k) + thickness(k) ) / 2.0 ! 2nd order
- end do
+ else
+ print *,'layerThicknessEdge not calculated!'
+ endif
-
-
! Calculate masks
-
call land_ice_calculate_mask(mesh, state, err)
@@ -538,7 +541,7 @@
!> results.
!
!-----------------------------------------------------------------------
- subroutine land_ice_tend_h_layers_fo_upwind(grid, normalVelocity, layerThicknessEdge, tend, dt, err)!{{{
+ subroutine land_ice_tend_h_layers_fo_upwind(grid, normalVelocity, layerThickness, tend, dt, err)!{{{
!-----------------------------------------------------------------
!
@@ -550,7 +553,7 @@
normalVelocity !< Input: velocity
real (kind=RKIND), dimension(:,:), intent(in) :: &
- layerThicknessEdge !< Input: thickness of each layer at edge
+ layerThickness !< Input: thickness of each layer
type (mesh_type), intent(in) :: &
grid !< Input: grid information
@@ -580,11 +583,11 @@
!
!-----------------------------------------------------------------
- integer :: iEdge, nEdges, cell1, cell2, nVertLevels, k
+ integer :: iEdge, nEdges, cell1, cell2, nVertLevels, k, CellDownwind, CellUpwind
integer, dimension(:,:), pointer :: cellsOnEdge
- real (kind=RKIND) :: flux
+ real (kind=RKIND) :: flux, VelSign
real (kind=RKIND), dimension(:), pointer :: dvEdge, areaCell, dcEdge
err = 0
@@ -600,14 +603,24 @@
do iEdge=1,nEdges
cell1 = cellsOnEdge(1,iEdge)
cell2 = cellsOnEdge(2,iEdge)
- do k=1, nVertLevels
- if ( (abs(normalVelocity(k, iEdge)) * dt/SecondsInYear) .gt. (0.5 * dcEdge(iEdge))) then
- print *, 'CFL violation at level, edge', k, iEdge
- endif
- flux = normalVelocity(k,iEdge) * dvEdge(iEdge) * layerThicknessEdge(k,iEdge)
- tend(k,cell1) = tend(k,cell1) - flux / areaCell(cell1)
- tend(k,cell2) = tend(k,cell2) + flux / areaCell(cell2)
- 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 (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)
+ do k=1, nVertLevels
+ if ( (abs(normalVelocity(k, iEdge)) * dt/SecondsInYear) .gt. (0.5 * dcEdge(iEdge))) then
+ print *, 'CFL violation at level, edge', k, iEdge
+ endif
+ flux = normalVelocity(k,iEdge) * dvEdge(iEdge) * layerThickness(k,cellUpwind)
+ tend(k,cell1) = tend(k,cell1) - flux / areaCell(cell1)
+ tend(k,cell2) = tend(k,cell2) + flux / areaCell(cell2)
+ end do
+ 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-05-21 20:07:29 UTC (rev 1924)
+++ branches/land_ice_projects/implement_core/src/core_land_ice/mpas_land_ice_time_integration_forwardeuler.F 2012-05-21 21:14:18 UTC (rev 1925)
@@ -45,7 +45,7 @@
type (state_type), pointer :: stateOld, stateNew
type (mesh_type), pointer :: mesh
- real (kind=RKIND), dimension(:), pointer :: thicknessOld, thicknessNew, thickness_tend, layerThicknessFractions
+ real (kind=RKIND), dimension(:), pointer :: thicknessOld, thicknessNew, layerThicknessFractions
real (kind=RKIND), dimension(:,:), pointer :: normalVelocityOld, normalVelocityNew, layerThicknessOld, layerThicknessNew, layerThickness_tend
real (kind=RKIND), dimension(:,:,:), pointer :: tracer_tendency, tracersNew, tracersOld
@@ -90,7 +90,6 @@
! Tendencies
layerThickness_tend => block % tend % layerThickness % array
- thickness_tend => block % tend % thickness % array
tracer_tendency => block % tend % tracers % array
</font>
</pre>