<p><b>mhoffman@lanl.gov</b> 2012-05-11 14:38:33 -0600 (Fri, 11 May 2012)</p><p>BRANCH COMMIT<br>
<br>
Created a mask module that calculates a mask for cells, vertices, and edges.<br>
Masks are integer bitmasks with bits for ice/no-ice, thin-ice/thick-ice, grounded/floating, etc. These are likely to evolve over time as the masks are used.<br>
Mask details are in the mpas_land_ice_mask.inc file, which needs to be included in any module that makes use of macros written to query the masks.<br>
Applied the edgeMask in the SIA velocity module.<br>
</p><hr noshade><pre><font color="gray">Modified: branches/land_ice_projects/implement_core/Makefile
===================================================================
--- branches/land_ice_projects/implement_core/Makefile 2012-05-11 15:55:34 UTC (rev 1898)
+++ branches/land_ice_projects/implement_core/Makefile 2012-05-11 20:38:33 UTC (rev 1899)
@@ -207,7 +207,8 @@
####################################################
RM = rm -f
-CPP = cpp -C -P -traditional
+# -C flag needed to be removed from CPP call to allow comments in the mask include file.
+CPP = cpp -P -traditional
RANLIB = ranlib
#########################
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-11 15:55:34 UTC (rev 1898)
+++ branches/land_ice_projects/implement_core/src/core_land_ice/Makefile 2012-05-11 20:38:33 UTC (rev 1899)
@@ -11,6 +11,7 @@
mpas_land_ice_tendency.o \
mpas_land_ice_time_integration.o \
mpas_land_ice_time_integration_forwardeuler.o \
+ mpas_land_ice_mask.o \
mpas_land_ice_vel.o \
mpas_land_ice_lifev.o \
mpas_land_ice_sia.o \
@@ -25,9 +26,11 @@
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_tendency.o: mpas_land_ice_mask.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-11 15:55:34 UTC (rev 1898)
+++ branches/land_ice_projects/implement_core/src/core_land_ice/Registry 2012-05-11 20:38:33 UTC (rev 1899)
@@ -21,6 +21,10 @@
namelist character land_ice_model config_advection FO-Upwind
% valid advection is FO-Upwind, None (currently only applies to thickness)
namelist real land_ice_model config_ice_density 900.0
+namelist real land_ice_model config_ocean_density 1000.0
+namelist real land_ice_model config_sealevel 0.0
+namelist real land_ice_model config_dynamic_thickness 100.0
+% Define the ice thickness below which dynamics are not calculated.
% 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
@@ -168,6 +172,7 @@
% Prognostic variables: read from input, saved in restart, and written to output
var persistent real thickness ( nCells Time ) 2 iro thickness state - -
var persistent real bedTopography ( nCells Time ) 2 iro bedTopography state - -
+% bedTopography is really a mesh variable now, but would be a state variable if isostasy is included.
var persistent real temperature ( nVertLevels nCells Time ) 2 iro temperature state tracers dynamics
var persistent real normalVelocity ( nVertLevels nEdges Time ) 2 iro normalVelocity state - -
Added: 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 (rev 0)
+++ branches/land_ice_projects/implement_core/src/core_land_ice/mpas_land_ice_mask.F 2012-05-11 20:38:33 UTC (rev 1899)
@@ -0,0 +1,249 @@
+!|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
+!
+! land_ice_mask
+!
+!> \MPAS land-ice mask calculations
+!> \author Matt Hoffman
+!> \date 10 May 2012
+!> \version SVN:$Id:$
+!> \details
+!> This module contains the routines for calculating masks for land ice
+!>
+!
+!-----------------------------------------------------------------------
+
+#include "mpas_land_ice_mask.inc"
+
+module land_ice_mask
+
+ use mpas_grid_types
+ use mpas_configure
+ use mpas_dmpar
+
+ implicit none
+ private
+ save
+
+ !--------------------------------------------------------------------
+ !
+ ! Public parameters
+ !
+ !--------------------------------------------------------------------
+
+ !--------------------------------------------------------------------
+ !
+ ! Public member functions
+ !
+ !--------------------------------------------------------------------
+ public :: land_ice_calculate_mask
+
+ !--------------------------------------------------------------------
+ !
+ ! Private module variables
+ !
+ !--------------------------------------------------------------------
+
+
+
+!***********************************************************************
+
+contains
+
+!***********************************************************************
+!
+! routine land_ice_calculate_mask
+!
+!> \brief Calculates masks for land ice
+!> \author Matt Hoffman
+!> \date 10 May 2012
+!> \version SVN:$Id$
+!> \details
+!> This routine Calculates masks for land ice.
+!
+!-----------------------------------------------------------------------
+
+ subroutine land_ice_calculate_mask(mesh, state, err)
+
+ !-----------------------------------------------------------------
+ !
+ ! input variables
+ !
+ !-----------------------------------------------------------------
+
+ type (mesh_type), intent(in) :: &
+ mesh !< Input: mesh information
+
+ !-----------------------------------------------------------------
+ !
+ ! input/output variables
+ !
+ !-----------------------------------------------------------------
+
+ type (state_type), intent(inout) :: &
+ state !< Input: state information
+
+ !-----------------------------------------------------------------
+ !
+ ! output variables
+ !
+ !-----------------------------------------------------------------
+
+ integer, intent(out) :: err !< Output: error flag
+
+ !-----------------------------------------------------------------
+ !
+ ! local variables
+ !
+ !-----------------------------------------------------------------
+
+
+ integer :: nCells, nVertices, nEdges
+
+ real(KIND=RKIND), dimension(:), pointer :: thickness, bedTopography
+
+ integer, dimension(:), pointer :: nEdgesOnCell, cellMask, vertexMask, edgeMask
+
+ integer, dimension(:,:), pointer :: cellsOnCell, cellsOnVertex, cellsOnEdge
+
+ integer :: i, j, iCell
+
+ logical :: isMargin
+
+ logical :: aCellOnVertexHasIce, aCellOnVertexHasThickIce, aCellOnVertexHasNoIce, aCellOnVertexIsFloating
+
+ logical :: aCellOnEdgeHasIce, aCellOnEdgeHasThickIce, aCellOnEdgeHasNoIce, aCellOnEdgeIsFloating
+
+
+ err = 0
+
+
+ ! Assign pointers and variables
+ cellMask => state % cellMask % array
+ thickness => state % thickness % array
+ bedTopography => state % bedTopography % array
+ vertexMask => state % vertexMask % array
+ edgeMask => state % edgeMask % array
+ nEdgesOnCell => mesh % nEdgesOnCell % array
+ cellsOnCell => mesh % cellsOnCell % array
+ cellsOnVertex => mesh % cellsOnVertex % array
+ cellsOnEdge => mesh % cellsOnEdge % array
+
+ nCells = mesh % nCells
+ nVertices = mesh % nVertices
+ nEdges = mesh % nEdges
+
+
+ ! Calculate cellMask values===========================
+ cellMask = 0
+
+ ! Bit 0: Identify points with ice
+ where (thickness > 0)
+ cellMask = ior(cellMask, MASK_VALUE_ICE)
+ end where
+
+ ! Bit 1: Identify points where the ice is above the ice dynamics limit
+ where ( thickness > config_dynamic_thickness )
+ cellMask = ior(cellMask, MASK_VALUE_THICK_ICE)
+ end where
+
+ ! 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.
+ where ( MASK_IS_ICE(cellMask) .and. (config_ice_density / config_ocean_density * thickness) .gt. (config_sealevel - bedTopography) )
+ cellMask = ior(cellMask, MASK_VALUE_GROUNDED)
+ end where
+ where ( MASK_IS_ICE(cellMask) .and. ( .not. MASK_IS_GROUNDED(cellMask) ) )
+ cellMask = ior(cellMask, MASK_VALUE_FLOATING)
+ end where
+
+ ! Bit 5: Identify the margin
+ ! For a cell, we define the margin as the first cell with ice (at least one neighbor is a non-ice cell)
+ do i=1,nCells
+ if MASK_IS_ICE(cellMask(i)) then
+ isMargin = .false.
+ do j=1,nEdgesOnCell(i) ! Check if any neighbors are non-ice
+ isMargin = ( isMargin .or. MASK_IS_NOT_ICE(cellMask(cellsOnCell(j,i))) )
+ enddo
+ 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!!!!!
+
+
+
+
+ ! 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
+ ! Bit 2: Grounded vertices have all neighboring cells with ice grounded (mutually exclusive with floating)
+ ! Bit 3: Floating vertices have at least one neighboring cell floating
+ ! Bit 4: G.L. vertices have at least one neighboring cell floating and at least one neighboring cell grounded (not yet implemented)
+ ! Bit 5: Vertices on margin are vertices with at least one neighboring cell with ice and at least one neighboring cell without ice
+ vertexMask = 0
+ do i = 1,nVertices
+ aCellOnVertexHasIce = .false.
+ aCellOnVertexHasThickIce = .false.
+ aCellOnVertexHasNoIce = .false.
+ aCellOnVertexIsFloating = .false.
+ do j = 1, mesh % vertexDegree ! vertexDegree is usually 3 (e.g. CVT mesh) but could be something else (e.g. 4 for quad mesh)
+ iCell = cellsOnVertex(j,i)
+ aCellOnVertexHasIce = (aCellOnVertexHasIce .or. MASK_IS_ICE(cellMask(iCell)))
+ aCellOnVertexHasThickIce = (aCellOnVertexHasThickIce .or. MASK_IS_THICK_ICE(cellMask(iCell)))
+ aCellOnVertexHasNoIce = (aCellOnVertexHasNoIce .or. MASK_IS_NOT_ICE(cellMask(iCell)))
+ aCellOnVertexIsFloating = (aCellOnVertexIsFloating .or. MASK_IS_FLOATING(cellMask(iCell)))
+ end do
+ vertexMask(i) = ior(vertexMask(i), MACRO_LOG_TO_INT(aCellOnVertexHasIce) * MASK_VALUE_ICE)
+ vertexMask(i) = ior(vertexMask(i), MACRO_LOG_TO_INT(aCellOnVertexHasThickIce) * MASK_VALUE_THICK_ICE)
+ vertexMask(i) = ior(vertexMask(i), MACRO_LOG_TO_INT(aCellOnVertexHasIce .and. aCellOnVertexHasNoIce) * MASK_VALUE_MARGIN) ! vertex with both 1+ ice cell and 1+ non-ice cell as neighbors
+ vertexMask(i) = ior(vertexMask(i), MACRO_LOG_TO_INT(aCellOnVertexIsFloating) * MASK_VALUE_FLOATING)
+ vertexMask(i) = ior(vertexMask(i), MACRO_LOG_TO_INT( MASK_IS_ICE(vertexMask(i)) .and. (.not. aCellOnVertexIsFloating)) * MASK_VALUE_GROUNDED) ! vertex that doesn't have any floating cells adjacent but does have ice
+ end do
+
+
+
+
+ ! Calculate edgeMask values based on cellMask values===========================
+ ! Bit 0: Edges with ice are ones with at least one adjacent cell with ice
+ ! Bit 1: Edges with thick ice are ones with at least one adjacent cell with thick ice
+ ! Bit 2: Grounded Edges have all neighboring cells with ice grounded (mutually exclusive with floating)
+ ! Bit 3: Floating Edges have at least one neighboring cell floating
+ ! Bit 4: G.L. Edges have one neighboring cell floating and one neighboring cell grounded (not yet implemented)
+ ! Bit 5: Edges on margin are vertices with one neighboring cell with ice and one neighboring cell without ice
+ edgeMask = 0
+ do i = 1,nEdges
+ aCellOnEdgeHasIce = .false.
+ aCellOnEdgeHasThickIce = .false.
+ aCellOnEdgeHasNoIce = .false.
+ aCellOnEdgeIsFloating = .false.
+ do j = 1, 2
+ iCell = cellsOnEdge(j,i)
+ aCellOnEdgeHasIce = (aCellOnEdgeHasIce .or. MASK_IS_ICE(cellMask(iCell)))
+ aCellOnEdgeHasThickIce = (aCellOnEdgeHasThickIce .or. MASK_IS_THICK_ICE(cellMask(iCell)))
+ aCellOnEdgeHasNoIce = (aCellOnEdgeHasNoIce .or. MASK_IS_NOT_ICE(cellMask(iCell)))
+ aCellOnEdgeIsFloating = (aCellOnEdgeIsFloating .or. MASK_IS_FLOATING(cellMask(iCell)))
+ end do
+ edgeMask(i) = ior(edgeMask(i), MACRO_LOG_TO_INT(aCellOnEdgeHasIce) * MASK_VALUE_ICE)
+ edgeMask(i) = ior(edgeMask(i), MACRO_LOG_TO_INT(aCellOnEdgeHasThickIce) * MASK_VALUE_THICK_ICE)
+ edgeMask(i) = ior(edgeMask(i), MACRO_LOG_TO_INT(aCellOnEdgeHasIce .and. aCellOnEdgeHasNoIce) * MASK_VALUE_MARGIN)
+ edgeMask(i) = ior(edgeMask(i), MACRO_LOG_TO_INT(aCellOnEdgeIsFloating) * MASK_VALUE_FLOATING)
+ edgeMask(i) = ior(edgeMask(i), MACRO_LOG_TO_INT( MASK_IS_ICE(edgeMask(i)) .and. (.not. aCellOnEdgeIsFloating)) * MASK_VALUE_GROUNDED)
+ end do
+
+
+ !--------------------------------------------------------------------
+
+ end subroutine land_ice_calculate_mask
+
+
+!***********************************************************************
+! Private subroutines:
+
+
+
+!***********************************************************************
+
+end module land_ice_mask
+
+!|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
+
Added: 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 (rev 0)
+++ branches/land_ice_projects/implement_core/src/core_land_ice/mpas_land_ice_mask.inc 2012-05-11 20:38:33 UTC (rev 1899)
@@ -0,0 +1,117 @@
+/*
+ Adapted from CISM's glide_mask.inc. Values not currently used have been commented with C-style comments using slash star
+ All mask parameters start with MASK_VALUE_ and all mask macros start with MASK_IS_
+
+ To use any of these definitions in a module, add #include "mpas_land_ice_mask.inc" to the very beginning (before the module declaration)
+*/
+
+
+
+/* Define Mask Values ====================================== */
+
+/*Bits 1:0 - Ice presence (0 if no ice) */
+#define MASK_VALUE_ICE 1
+/*NOTE: If bit 2 is activated, bit 1 must be activated*/
+#define MASK_VALUE_THICK_ICE 2
+/* #define GLIDE_ICE_PRESENCE_BITS 3 */
+
+/*Bits 4:2 - Type of base (Land or ocean - grounding line has both bits on)
+ 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.
+ The 16 bit specifies grounding line - CISM sets this up so that those points are treated as grounded as well. */
+#define MASK_VALUE_GROUNDED 4
+#define MASK_VALUE_FLOATING 8
+#define MASK_VALUE_GROUNDING_LINE 16
+/* #define GLIDE_BASE_TYPE_BITS 28 */
+
+/*Bit 5: Identifies a margin (a cell with nonzero thickness adjacent on any side to a cell with zero thickness)
+ 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: Identifies a dirichlet condition
+ !The velocity at points marked with this bit should be held constant.
+ #define GLIDE_MASK_DIRICHLET_BC 64
+
+ !Bit 7: Identifies a computational domain boundary. These are normally just activated on the edges of the domain,
+ !unless there is a domain decomposition (in which case they may be missing)
+ #define GLIDE_MASK_COMP_DOMAIN_BND 128
+*/
+
+
+
+
+/* Define Mask Macros ============================================= */
+
+/* NOTE: All Macros will return LOGICAL types. If you want the integer form (0/1), use the transfer function to 'typecast' from logical to integer (see macro below). Transfer will always convert false to 0, but may convert true to any nonzero value. However, values other than -1 or 1 have not been encountered. */
+
+/* A simple macro to convert from logical to integer */
+#define MACRO_LOG_TO_INT(logicalvalue) (abs(transfer(logicalvalue, 1)))
+
+/* A simple macro to do integer 1/0 NOT operations */
+#define MACRO_INOT(intvalue) (abs(intvalue-1))
+
+
+
+/* ===Checks for an iceless square: */
+/*Check for a lack of ice:*/
+#define MASK_IS_NOT_ICE(mask) (iand(mask, MASK_VALUE_ICE) == 0)
+
+
+/* ===Checks for ice: */
+/*Checks for the presence of any ice, dynamic or not*/
+#define MASK_IS_ICE(mask) (iand(mask, MASK_VALUE_ICE) == MASK_VALUE_ICE)
+
+/*Checks for the presence of ice that is above the ice dynamics limit*/
+#define MASK_IS_THICK_ICE(mask) (iand(mask, MASK_VALUE_THICK_ICE) == MASK_VALUE_THICK_ICE)
+
+/*Checks for the presence of ice that is below the ice dynamics limit*/
+#define MASK_IS_THIN_ICE(mask) (iand(mask, MASK_VALUE_THICK_ICE) == 0 .and. MASK_IS_ICE(mask))
+
+/*Checks for any ice, dynamic or not, that is on an ice shelf.*/
+#define MASK_IS_FLOATING(mask) (iand(mask, MASK_VALUE_FLOATING) == MASK_VALUE_FLOATING .and. MASK_IS_ICE(mask))
+
+/*Checks for any ice, dynamic or not, that is grounded*/
+#define MASK_IS_GROUNDED(mask) (iand(mask, MASK_VALUE_GROUNDED) == MASK_VALUE_GROUNDED .and. MASK_IS_ICE(mask))
+
+/*Checks whether this is an ice margin (first cell on glacier side from non-ice cell)*/
+#define MASK_IS_MARGIN(mask) (iand(mask, MASK_VALUE_MARGIN) == MASK_VALUE_MARGIN)
+
+
+
+
+/* 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)
+
+!Check for land with no ice:
+#define MASK_IS_LAND(mask) abs(transfer(iand(mask, MASK_VALUE_LAND_BASED) == MASK_VALUE_LAND_BASED .and. MASK_IS_NOT_ICE(mask)), 1)
+
+!Checks for any ice, dynamic or not, that is on the grounding line
+ #define MASK_IS_GROUNDING_LINE(mask) (iand(mask, 17) == 17)
+
+ Skipping these for now until we resolve how to deal with the G.L.
+ !Checks for any ice, dynamic or not, that is either floating *or* on the grounding line
+ #define GLIDE_IS_FLOAT_OR_GNDLINE(mask) (iand(mask, 24) > 0 .and. GLIDE_HAS_ICE(mask))
+
+ !Checks for any ice, dynamic or not, that is either grounded *or* on the grounding line
+ #define GLIDE_IS_GROUND_OR_GNDLINE(mask) (iand(mask, 20) > 0 .and. GLIDE_HAS_ICE(mask))
+
+!Checks whether this is a margin in contact with the ocean, floating ice or open ocean
+ #define MASK_IS_MARINE_MARGIN(mask) (MASK_IS_MARGIN(mask) .and. GLIDE_IS_FLOAT_OR_GNDLINE(mask))
+
+!Checks whether this is a calving margin in contact with the ocean
+!41 = 32+8+1 = margin + ocean-based + ice
+#define MASK_IS_CALVING(mask) (iand(mask, 41) == 41)
+
+!Checks whether this is a land margin
+!37 = 32 + 4 + 1 = margin + land-based + ice
+#define MASK_IS_LAND_MARGIN(mask) (iand(mask, 37) == 37)
+
+ !Checks whether a dirichlet boundary has been defined at this point
+ #define GLIDE_IS_DIRICHLET_BOUNDARY(mask) (iand(mask, GLIDE_MASK_DIRICHLET_BC) == GLIDE_MASK_DIRICHLET_BC)
+ !Checks whether we are at the edge of the computational domain *and* there is ice in this square
+ #define GLIDE_IS_COMP_DOMAIN_BND(mask) (iand(mask, 129) == 129)
+*/
+
+
+
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-11 15:55:34 UTC (rev 1898)
+++ branches/land_ice_projects/implement_core/src/core_land_ice/mpas_land_ice_sia.F 2012-05-11 20:38:33 UTC (rev 1899)
@@ -12,6 +12,8 @@
!
!-----------------------------------------------------------------------
+#include "mpas_land_ice_mask.inc"
+
module land_ice_sia
use mpas_grid_types
@@ -253,10 +255,10 @@
real (kind=RKIND), dimension(:), pointer :: dcEdge
real (kind=RKIND), dimension(:,:), pointer :: layerThickness, 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
- logical edgeIsIce
real rhoi !!!! THIS SHOULD BE A PHYSICAL PARAMETER ELSEWHERE
err = 0
@@ -273,6 +275,7 @@
normalVelocity => state % normalVelocity % array
thickness => state % thickness % array
layerThickness => state % layerThickness % array
+ edgeMask => state % edgeMask % array
rhoi = config_ice_density
@@ -289,12 +292,9 @@
! Loop over edges
do iEdge = 1, nEdges
- ! Determine if ice
- call is_edge_ice(thickness, cellsOnEdge, iEdge, edgeIsIce)
! Only calculate velocity for edges that are part of the ice sheet.
- ! This should be a mask calculated external this module.
! Also, the velocity calculation should be valid for non-ice edges (i.e. returns 0).
- if (edgeIsIce .eqv. .true.) then
+ 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)) ) &
@@ -319,35 +319,9 @@
! private subroutines
- subroutine is_edge_ice(thickness, cellsOnEdge, iEdge, edgeIsIce)
- ! Determine if an edge is part of the ice sheet.
- ! This should probably be calculated external to this module as a mask.
- ! (and could be more sophisticated than this)
+
- implicit none
- real (kind=RKIND), dimension(:), intent(in) :: thickness
- integer, dimension(:,:), intent(in) :: cellsOnEdge
- integer, intent(in) :: iEdge
- logical, intent(out) :: edgeIsIce
-
- ! Local variables
- integer :: cell1, cell2
-
- cell1 = cellsOnEdge(1,iEdge)
- cell2 = cellsOnEdge(2,iEdge)
- ! Check if an edge has at least one adjacent cell with ice
- ! (to turn off advancing margins with FO thickness advection, replace the .or. with .and.
- if ( (thickness(cell1).gt.0.0) .or. (thickness(cell2).gt.0.0) ) then
- edgeIsIce = .true.
- else
- edgeIsIce = .false.
- endif
-
- end subroutine is_edge_ice
-
-
-
!***********************************************************************
end module land_ice_sia
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-11 15:55:34 UTC (rev 1898)
+++ branches/land_ice_projects/implement_core/src/core_land_ice/mpas_land_ice_tendency.F 2012-05-11 20:38:33 UTC (rev 1899)
@@ -19,6 +19,7 @@
use mpas_configure
use mpas_constants
use mpas_dmpar
+ use land_ice_mask
implicit none
private
@@ -214,19 +215,9 @@
!\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.
- !\todo Calculate masks. Right now sia driver and Forward Euler time integration have need of a mask and both do it on their own. We need to define how the masks should work and implement them here.
- ! simple vertex mask needed by LifeV
- do i = 1,nVertices
- keep_vertex = 0
- do k = 1, 3
- iCell = cellsOnVertex(k,i)
- if (thickness(iCell) .gt. 0.01) then
- keep_vertex = 1
- endif
- end do
- vertexMask(i) = keep_vertex
- end do
+!!!!!!!!!!!!!!!!!! !\todo Calculate masks. Right now sia driver and Forward Euler time integration have need of a mask and both do it on their own. We need to define how the masks should work and implement them here.
+ call land_ice_calculate_mask(mesh, state, err)
end subroutine land_ice_diagnostic_solve
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-11 15:55:34 UTC (rev 1898)
+++ branches/land_ice_projects/implement_core/src/core_land_ice/mpas_land_ice_time_integration_forwardeuler.F 2012-05-11 20:38:33 UTC (rev 1899)
@@ -149,7 +149,19 @@
call mpas_timer_start("calc. diagnostic vars except vel")
call land_ice_diagnostic_solve(mesh, stateNew, err) ! perhaps velocity solve should move in here.
+ ! update halos on masks (margins need neighbor info, so the outermost halo levels could be wrong
+ call mpas_dmpar_exch_halo_field1d_integer(dminfo, stateNew % cellMask % array, &
+ mesh % nCells, &
+ block % parinfo % cellsToSend, block % parinfo % cellsToRecv)
+ call mpas_dmpar_exch_halo_field1d_integer(dminfo, stateNew % vertexMask % array, &
+ mesh % nVertices, &
+ block % parinfo % verticesToSend, block % parinfo % verticesToRecv)
+ call mpas_dmpar_exch_halo_field1d_integer(dminfo, stateNew % edgeMask % array, &
+ mesh % nEdges, &
+ block % parinfo % edgesToSend, block % parinfo % edgesToRecv)
+
! Determine if the vertex mask changed during this time step for this block (needed for LifeV)
+ ! \todo: there may be some aspects of the mask that are ok change for LifeV, but for now just check the whole thing.
if ( sum(stateNew % vertexMask % array - stateOld % vertexMask % array) .ne. 0 ) then
blockVertexMaskChanged = 1
else
</font>
</pre>