<p><b>mhoffman@lanl.gov</b> 2013-02-01 10:40:00 -0700 (Fri, 01 Feb 2013)</p><p>BRANCH COMMIT - land ice<br>
<br>
Implementing ice shelf advance using the Calving Front Boundary Condition approach described by Albrecht et al. 2012 (www.the-cryosphere.net/5/35/2011/). This prevents the formation of unrealistic long, thin ice shelves advancing dispersively which happens without a special treatment of ice shelf advance. Basically, the technique tracks partially filled cells along the ice shelf margin, preventing advance to the next cell until a critical thickness is filled. This is a necessary feature to get a realistic implementation of even simple calving laws.<br>
<br>
This option is turned on with Namelist option: config_marine_advance='CFBC'. To turn it off, use 'basic' (or anything else at the moment). 'basic' is the default for now.<br>
<br>
In my implementation, all marginal ice shelf cells are dynamically inactive (the same as 'thin ice'). The reference height for a partially filled cell is set to the mean height of all cells feeding it (in a FO Upwinding sense). (In the future this may be a flux-weighted mean.) Once the cell fills, any residual flux is divided evenly between open ocean cells in front of it. Velocities are then back-calculated on these newly active edges to be consistent with the flux moved across them (velocities are needed for tracer advection). This entire process is an additional step during thickness tendency calculation.<br>
<br>
Changes were also made to the mask implementation to be consistent with the CFBC conventions when CFBC is on, and bug fixes/improvements were necessary in mpas_init_block(), land_ice_sia_solve(), and land_ice_time_integrator_ForwardEuler().<br>
</p><hr noshade><pre><font color="gray">Modified: branches/land_ice_projects/implement_core/src/core_land_ice/Registry
===================================================================
--- branches/land_ice_projects/implement_core/src/core_land_ice/Registry 2013-02-01 06:19:48 UTC (rev 2409)
+++ branches/land_ice_projects/implement_core/src/core_land_ice/Registry 2013-02-01 17:40:00 UTC (rev 2410)
@@ -30,6 +30,8 @@
namelist character land_ice_model config_forcing_frequency Static
% Valid options are Static, Annual
namelist logical land_ice_model config_allow_additional_advance true
+namelist character land_ice_model config_marine_advance basic
+% valid options are CFBC and basic
% The following options are used by the ocean core and may be useful in the future
%namelist logical land_ice_model config_h_ScaleWithMesh false
Modified: branches/land_ice_projects/implement_core/src/core_land_ice/mpas_land_ice_mask.F
===================================================================
--- branches/land_ice_projects/implement_core/src/core_land_ice/mpas_land_ice_mask.F 2013-02-01 06:19:48 UTC (rev 2409)
+++ branches/land_ice_projects/implement_core/src/core_land_ice/mpas_land_ice_mask.F 2013-02-01 17:40:00 UTC (rev 2410)
@@ -221,6 +221,7 @@
where ( thickness > config_dynamic_thickness )
cellMask = ior(cellMask, MASK_VALUE_THICK_ICE)
end where
+ ! see below for an additional check
! Bit 2 & 3: Is it grounded? (ice thickness equal to floatation is considered floating)
! Note: we may want to make these the same bit, and say that all ice is either floating or grounded - but depends on how we define the grounding line.
@@ -242,13 +243,18 @@
cellMask(i) = ior(cellMask(i), MACRO_LOG_TO_INT(isMargin) * MASK_VALUE_MARGIN)
endif
enddo
+ ! TODO Probably also want to identify the margin between thin and thick ice!!!!!
+ ! Bit 1: redux
+ ! If we are using CFBC for marine advance, then we want to make the last ice shelf cell dynamically inactive
+ ! TODO Do we not want to do this at the initial time? If we do it then we lose dynamics on the last cell
+ if (config_marine_advance == 'CFBC') then
+ where ( (MASK_IS_FLOATING(cellMask)) .and. (MASK_IS_MARGIN(cellMask)) )
+ cellMask = iand(cellMask, not(MASK_VALUE_THICK_ICE)) ! Turn off the thick ice flag
+ end where
+ end if
- ! \todo Probably also want to identify the margin between thin and thick ice!!!!!
-
-
-
! Calculate vertexMask values based on cellMask values===========================
! Bit 0: Vertices with ice are ones with at least one adjacent cell with ice
! Bit 1: Vertices with thick ice are ones with at least one adjacent cell with thick ice
Modified: branches/land_ice_projects/implement_core/src/core_land_ice/mpas_land_ice_mask.inc
===================================================================
--- branches/land_ice_projects/implement_core/src/core_land_ice/mpas_land_ice_mask.inc 2013-02-01 06:19:48 UTC (rev 2409)
+++ branches/land_ice_projects/implement_core/src/core_land_ice/mpas_land_ice_mask.inc 2013-02-01 17:40:00 UTC (rev 2410)
@@ -27,7 +27,13 @@
The type of margin is determined by whether the ice here is grounded or floating (in conjunction with those flags) */
#define MASK_VALUE_MARGIN 32
+
+
/* ====Not used, currently:
+
+! Bit 6: if the bed is below sea level ! Decided this was not worth vaculating every time step since it doesn't change currently, and can easily be determined
+#define MASK_VALUE_OCEAN_BASED 64
+
!Bit 6: Identifies a dirichlet condition
!The velocity at points marked with this bit should be held constant.
#define GLIDE_MASK_DIRICHLET_BC 64
@@ -83,6 +89,16 @@
+
+/* I wrote these but decided they were not worth having around
+! Check for ocean-based bed
+#define MASK_IS_OCEAN_BASED(mask) (iand(mask, MASK_VALUE_OCEAN_BASED) == MASK_VALUE_OCEAN_BASED)
+
+! check for open ocean with no ice
+#define MASK_IS_OPEN_OCEAN(mask) (not(MASK_IS_ICE(mask)) .and. (MASK_IS_OCEAN_BASED(mask)) )
+*/
+
+
/* Ignoring these for now (need to be checked/modified before using!):
!Check for open ocean with no ice:
#define MASK_IS_OCEAN(mask) abs(transfer(iand(mask, MASK_VALUE_OCEAN_BASED) == MASK_VALUE_OCEAN_BASED .and. MASK_IS_NOT_ICE(mask)), 1)
Modified: branches/land_ice_projects/implement_core/src/core_land_ice/mpas_land_ice_mpas_core.F
===================================================================
--- branches/land_ice_projects/implement_core/src/core_land_ice/mpas_land_ice_mpas_core.F 2013-02-01 06:19:48 UTC (rev 2409)
+++ branches/land_ice_projects/implement_core/src/core_land_ice/mpas_land_ice_mpas_core.F 2013-02-01 17:40:00 UTC (rev 2410)
@@ -154,7 +154,7 @@
use mpas_rbf_interpolation
use mpas_vector_reconstruction
use land_ice_vel
- use land_ice_tendency, only: land_ice_diagnostic_solve
+ use land_ice_tendency, only: land_ice_diagnostic_solve, land_ice_diagnostic_solve_using_velocity
use land_ice_mask
use mpas_ocn_tracer_advection
use ocn_advection
@@ -213,8 +213,24 @@
! \todo: skip this step if a velocity field was supplied in the I.C. input file
call mpas_timer_start("initial velocity calculation")
call land_ice_vel_solve(mesh, state, err)
+
+ ! TODO Iniitalize state may need to move out of block_init and into plain ol' init because a halo update may be needed after the velo solve before the mpas_reconstruct call
+ !!!! 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 the initial time - do this after velocity halo update in case velocities on the 1-halo edge are wrong (depends on velocity solver)
+ call mpas_reconstruct(mesh, state % normalVelocity % array,&
+ state % uReconstructX % array, &
+ state % uReconstructY % array, &
+ state % uReconstructZ % array, &
+ state % uReconstructZonal % array, &
+ state % uReconstructMeridional % array &
+ )
call mpas_timer_stop("initial velocity calculation")
+ call land_ice_diagnostic_solve_using_velocity(mesh, state, err) ! Some diagnostic variables require velocity to compute
if(err == 1) then
print *, "An error has occurred in mpas_init_block. Aborting..."
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 2013-02-01 06:19:48 UTC (rev 2409)
+++ branches/land_ice_projects/implement_core/src/core_land_ice/mpas_land_ice_sia.F 2013-02-01 17:40:00 UTC (rev 2410)
@@ -296,9 +296,9 @@
! Loop over edges
do iEdge = 1, nEdges
- ! Only calculate velocity for edges that are part of the ice sheet.
+ ! Only calculate velocity for edges that are part of the dynamic ice sheet.(thick ice)
! Also, the velocity calculation should be valid for non-ice edges (i.e. returns 0).
- if ( MASK_IS_ICE(edgeMask(iEdge)) ) then
+ if ( MASK_IS_THICK_ICE(edgeMask(iEdge)) ) then
cell1 = cellsOnEdge(1,iEdge)
cell2 = cellsOnEdge(2,iEdge)
! Calculate slope at edge
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-01 06:19:48 UTC (rev 2409)
+++ branches/land_ice_projects/implement_core/src/core_land_ice/mpas_land_ice_tendency.F 2013-02-01 17:40:00 UTC (rev 2410)
@@ -41,7 +41,8 @@
!--------------------------------------------------------------------
public :: land_ice_tendency_h, &
land_ice_tendency_tracers, &
- land_ice_diagnostic_solve
+ land_ice_diagnostic_solve, &
+ land_ice_diagnostic_solve_using_velocity
!--------------------------------------------------------------------
!
@@ -80,8 +81,9 @@
type (mesh_type), intent(in) :: &
mesh !< Input: grid information
- type (state_type), intent(in) :: &
- state !< Input: state to use to calculate tendency
+ type (state_type), intent(inout) :: &
+ state !< Input: state to use to calculate tendency (old time level)
+ ! Note: state needs to be inout (rather than just in) so that adjust_marine_boundary_fluxes can modify it.
real (kind=RKIND), intent(in) :: &
dt !< Input: dt
@@ -128,8 +130,12 @@
! state % layerThickness % array, state % thickness % array, thickness_tend, dt, err)
call land_ice_tend_h_layers_fo_upwind(mesh, state % normalVelocity % array, &
- state % layerThickness % array, layerThickness_tend, dt, allowableDt, err)
+ state % layerThickness % array, state % edgeMask % array, layerThickness_tend, dt, allowableDt, err)
+ ! Now adjust for marine boundary advance if needed
+ call adjust_marine_boundary_fluxes(mesh, state % thickness % array, state % layerThicknessEdge % array, &
+ state % cellMask % array, state % bedTopography % array, dt, layerThickness_tend, state % normalVelocity % array, err)
+
!case ('FCT') !===================================================
! MJH TEMP TRYING IT WITH THICKNESS
! stateOld % sup_thickness % array(1,:,:) = stateOld % layerThickness % array !!!!! MJH TEMP
@@ -649,7 +655,7 @@
!> results.
!
!-----------------------------------------------------------------------
- subroutine land_ice_tend_h_layers_fo_upwind(grid, normalVelocity, layerThickness, tend, dt, MinOfMaxAllowableDt, err)!{{{
+ subroutine land_ice_tend_h_layers_fo_upwind(grid, normalVelocity, layerThickness, edgeMask, tend, dt, MinOfMaxAllowableDt, err)!{{{
!-----------------------------------------------------------------
!
@@ -657,14 +663,17 @@
!
!-----------------------------------------------------------------
+ type (mesh_type), intent(in) :: &
+ grid !< Input: grid information
+
real (kind=RKIND), dimension(:,:), intent(in) :: &
normalVelocity !< Input: velocity
real (kind=RKIND), dimension(:,:), intent(in) :: &
layerThickness !< Input: thickness of each layer
- type (mesh_type), intent(in) :: &
- grid !< Input: grid information
+ integer, dimension(:), intent(in) :: &
+ edgeMask !< Input: mask on edges
real (kind=RKIND) :: dt
@@ -713,9 +722,10 @@
do iEdge=1,nEdges
cell1 = cellsOnEdge(1,iEdge)
cell2 = cellsOnEdge(2,iEdge)
- if ( (layerThickness(1,cell1) .gt. 0.0) .or. (layerThickness(1,cell2) .gt. 0.0) ) then
+! 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
+ if ( MASK_IS_THICK_ICE(edgeMask(iEdge)) ) then
VelSignSum = 0.0
do k=1, nVertLevels
@@ -754,10 +764,235 @@
print *, ' Maximum allowable time step (yrs) on this processor is ~', MinOfMaxAllowableDt
+
!--------------------------------------------------------------------
end subroutine land_ice_tend_h_layers_fo_upwind!}}}
+
+
+
+!***********************************************************************
+!
+! subroutine adjust_marine_boundary_fluxes
+!
+!> \brief Adjusts thickness tendency for marine margin (ice shelf front) cells
+!> \author Matt Hoffman
+!> \date 01 February 2013
+!> \version SVN:$Id$
+!> \details
+!> Adjusts thickness tendency after performing a special treatment of marine
+!> margin (ice shelf front) cells to ensure a Calving Flux Boundary Condition (CFBC).
+!> This is based on the approach described by Albrecht et al. 2012 (www.the-cryosphere.net/5/35/2011/).
+!> This prevents the formation of unrealistic long, thin ice shelves advancing dispersively,
+!> which happens without a special treatment of ice shelf advance. Basically, the
+!> technique tracks partially filled cells along the ice shelf margin, preventing
+!> advance to the next cell until a critical thickness is filled.
+!> In my implementation, all marginal ice shelf cells are dynamically inactive
+!> (the same as 'thin ice'). The reference height for a partially filled cell
+!> is set to the mean height of all cells feeding it (in a FO Upwinding sense).
+!> (In the future this may be a flux-weighted mean.) Once the cell fills, any
+!> residual flux is divided evenly between open ocean cells in front of it.
+!> Velocities are then back-calculated on these newly active edges to be consistent
+!> with the flux moved across them (these velocities are needed for tracer advection).
+!> This option is turned on with Namelist option: config_marine_advance='CFBC'.
+!> To turn it off, use 'basic' (or anything else at the moment).
+!> Note: there also is CFBC-specific code in mpas_land_ice_mask.F, subroutine land_ice_calculate_mask.
+
+!
+!-----------------------------------------------------------------------
+ subroutine adjust_marine_boundary_fluxes(mesh, thickness, layerThicknessEdge, cellMask, bedTopography, deltat, layerThickness_tend, normalVelocity, err)!{{{
+
+ !-----------------------------------------------------------------
+ !
+ ! input variables
+ !
+ !-----------------------------------------------------------------
+
+ type (mesh_type), intent(in) :: &
+ mesh !< Input: grid information
+
+ real (kind=RKIND), dimension(:), intent(in) :: &
+ thickness !< Input: thickness of each cell
+
+ real (kind=RKIND), dimension(:,:), intent(in) :: &
+ layerThicknessEdge !< Input: thickness of each layer on edges
+
+ integer, dimension(:), intent(in) :: &
+ cellMask !< Input: mask on cells
+
+ real (kind=RKIND), dimension(:), intent(in) :: &
+ bedTopography !< Input: thickness of each cell
+
+ real (kind=RKIND), intent(in) :: deltat !< Input: time step
+
+ !-----------------------------------------------------------------
+ !
+ ! input/output variables
+ !
+ !-----------------------------------------------------------------
+
+ real (kind=RKIND), dimension(:,:), intent(inout) :: &
+ layerThickness_tend !< Input: layer thickness tendency; modified on advancing marine margins
+
+ real (kind=RKIND), dimension(:,:), intent(inout) :: &
+ normalVelocity !< Input: velocity of edges; modified on advancing marine margins
+ !< Note: the modified velocity will not be written out because the output of its time level already occurred. However the modified velocity is needed for tracer advection.
+
+ !-----------------------------------------------------------------
+ !
+ ! output variables
+ !
+ !-----------------------------------------------------------------
+
+ integer, intent(out) :: err !< Output: error flag
+
+ !-----------------------------------------------------------------
+ !
+ ! local variables
+ !
+ !-----------------------------------------------------------------
+
+ ! pointers to mesh arrays
+ integer, dimension(:), pointer :: nEdgesOnCell
+ integer, dimension(:,:), pointer :: cellsOnCell, edgesOnCell, cellsOnEdge
+ real (kind=RKIND), dimension(:), pointer :: areaCell, dvEdge
+ ! counters, mesh variables, index variables
+ integer :: iCell, nCells, iEdge, neighborIndex, edgeIndex, k, nVertLevels, jEdge, kEdge
+ ! stuff for making calculations
+ real (kind=RKIND) :: numUpwindNeighbors, sumUpwindThickness, numOceanNeighbors, Href, projectedThickness, ExtraVolume, tendencyLayerAdjustment, VelSign, minHref, maxHref
+
+
+ ! Only execute the contents of this subroutine if CFBC is turned on.
+ if (config_marine_advance == 'CFBC') then
+
+ !write(0,*), 'Doing CFBC'
+
+ ! === Initialize variables and pointers ===
+
+ nCells = mesh % nCells
+ nVertLevels = mesh % nVertLevels
+
+ nEdgesOnCell => mesh % nEdgesOnCell % array
+ cellsOnCell => mesh % cellsOnCell % array
+ areaCell => mesh % areaCell % array
+ dvEdge => mesh % dvEdge % array
+ edgesOnCell => mesh % edgesOnCell % array
+ cellsOnEdge => mesh % cellsOnEdge % array
+
+ minHref = 99999.0
+ maxHref = -1.0
+
+ ! === Loop over cells to find marine margin cells ===
+ ! (I am assuming all marine margin cells require CFBC treatment -
+ ! - the occurrence of an exactly filled marine margin cell is not considered.
+ ! Thus the last floating cell along the margin will never be dynamically active with this scheme.)
+ do iCell = 1, nCells ! Need to solve over entire block to get correct values on the 0-halo
+ ! Identify where the partially filled marine margin cells are that will need to be dealt with
+ if ( (MASK_IS_FLOATING(cellMask(iCell))) .and. (MASK_IS_MARGIN(cellMask(iCell))) ) then
+
+ ! === Calculate the reference height from the edges feeding this cell ===
+ numUpwindNeighbors = 0.0
+ sumUpwindThickness = 0.0
+ numOceanNeighbors = 0.0
+ do iEdge = 1, nEdgesonCell(iCell)
+ edgeIndex = edgesOnCell(iEdge, iCell)
+ ! check which edges have velocity into this cell ! + velocity goes from index 1 to 2 in the cellsOnEdge array.
+ ! Note: I have left a commented out version that does a more robust check for net flux over all layers into the cell
+! if ( ((sum(normalVelocity(:,edgeIndex)*layerThicknessEdge(:, edgeIndex)) > 0.0) .and. (cellsOnEdge(2, edgeIndex) == iCell)) &
+! .or. ((sum(normalVelocity(:,edgeIndex)*layerThicknessEdge(:, edgeIndex)) < 0.0) .and. (cellsOnEdge(1, edgeIndex) == iCell)) ) then
+ ! Just check the surface velocity for directionality
+ if ( ((normalVelocity(1,edgeIndex) > 0.0) .and. (cellsOnEdge(2, edgeIndex) == iCell)) &
+ .or. ((normalVelocity(1,edgeIndex) < 0.0) .and. (cellsOnEdge(1, edgeIndex) == iCell)) ) then
+ numUpwindNeighbors = numUpwindNeighbors + 1.0
+ sumUpwindThickness = sumUpwindThickness + sum(layerThicknessEdge(:, edgeIndex))
+ end if
+
+ ! As long as we have an edge loop here, count up the number of open-ocean neighbors,
+ ! which will be needed if we have to move any residual mass forward to open ocean neighbors.
+ neighborIndex = cellsOnCell(iEdge, iCell)
+ if ( MASK_IS_NOT_ICE(cellMask(neighborIndex)) .and. (bedTopography(neighborIndex) < config_sealevel) ) then
+ numOceanNeighbors = numOceanNeighbors + 1.0
+ endif
+ end do
+ if (numUpwindNeighbors == 0.0) then
+ write(0,*) 'Error: in CFBC calculation an ice shelf margin cell has no ice neighbors!'
+ err = 1
+ else
+ !write(0,*) 'ok: margin cell has x neighbors', numUpwindNeighbors
+ endif
+ Href = sumUpwindThickness / numUpwindNeighbors
+ !write(0,*) 'Href=', Href
+ ! TODO make this a flux averaged reference height
+ minHref = min(minHref, Href)
+ maxHref = max(maxHref, Href)
+
+ ! === Calculate if/how much this cell will be overfilled by the new tendency ===
+ projectedThickness = thickness(iCell) + sum(layerThickness_tend(:,iCell)) * (deltat / SecondsInYear) ! The new thickness the cell will have after time integration
+ if (projectedThickness > Href) then
+ ! === Calculate residual thickness and move forward into open ocean cells ===
+ ExtraVolume = (projectedThickness - Href) * areaCell(iCell) ! units = m^3
+
+ do iEdge = 1, nEdgesonCell(iCell)
+ neighborIndex = cellsOnCell(iEdge, iCell)
+ if ( MASK_IS_NOT_ICE(cellMask(neighborIndex)) .and. (bedTopography(neighborIndex) < config_sealevel) ) then ! find open ocean neighbors
+
+ if (numOceanNeighbors<1) then
+ ! TODO Deal with this situation without a fatal error.
+ write(0,*) 'Error: in CFBC calculation an ice shelf margin cell needs to move residual mass forward but has no open-ice neighbors!'
+ err = 1
+ end if
+
+ ! Pass extra mass to neighboring ice-free (ocean only) cells by adjusting tendency in those cells
+ ! I am spreading the extra ice volume evenly amongst all vertical layers. This is a good approximation for SSA flow.
+ ! TODO this could be weighted by the edge length of the various neighbors
+ tendencyLayerAdjustment = ExtraVolume / numOceanNeighbors / areaCell(neighborIndex) / nVertLevels / (deltat / SecondsInYear)
+
+ ! Add the tendency to the neighbor. (Need to add it to the existing value
+ ! in case other overfilled margin cells have already contributed to this neighbor.)
+ layerThickness_tend(:,neighborIndex) = layerThickness_tend(:,neighborIndex) + tendencyLayerAdjustment
+ ! Lower the tendency in the overfilled cell by the same amount you are passing forward
+ layerThickness_tend(:,iCell) = layerThickness_tend(:,iCell) - tendencyLayerAdjustment
+
+ ! === Now calculate a consistent velocity on the edges passing mass to the new cells ===
+ ! Find the edge index (have to loop through edgesOnCell for both cells until we find the matching edge)
+ do jEdge = 1, nEdgesOnCell(neighborIndex)
+ do kEdge = 1, nEdgesOnCell(iCell)
+ if (edgesOnCell(jEdge, neighborIndex) == edgesOnCell(kEdge, iCell)) then
+ edgeIndex = edgesOnCell(kEdge, iCell)
+
+ if (cellsOnEdge(1, edgeIndex) == iCell) then
+ ! + velocity goes from index 1 to 2 in the cellsOnEdge array.
+ VelSign = 1.0
+ else
+ VelSign = -1.0
+ endif
+ do k = 1, nVertLevels
+ ! These velocities should have been 0 - we are now assigning a velocity
+ ! consistent with the flux we moved (using FO upwinding which is appropriate for advance)
+ normalVelocity(k,edgeIndex) = tendencyLayerAdjustment / &
+ (dvEdge(edgeIndex) * layerThicknessEdge(k, edgeIndex)) &
+ * areaCell(neighborIndex) * VelSign
+ end do
+ end if
+ end do
+ end do
+
+
+ end if ! open ocean neighbor
+ end do ! loop over edges on cell
+ end if ! overfilled cell
+ end if ! floating margin cell
+ end do ! cell loop
+
+ !write(0,*) 'CFBC: min, max href=', minHref, maxHref ! Useful for debugging.
+
+ end if ! to do CFBC or not
+
+ end subroutine adjust_marine_boundary_fluxes
+
+
+
end module land_ice_tendency
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 2013-02-01 06:19:48 UTC (rev 2409)
+++ branches/land_ice_projects/implement_core/src/core_land_ice/mpas_land_ice_time_integration_forwardeuler.F 2013-02-01 17:40:00 UTC (rev 2410)
@@ -295,7 +295,7 @@
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
+ call land_ice_diagnostic_solve_using_velocity(mesh, stateNew, err) ! Some diagnostic variables require velocity to compute
block => block % next
end do
call mpas_timer_stop("calc. diagnostic vars except vel")
</font>
</pre>