<p><b>mhoffman@lanl.gov</b> 2013-02-20 09:10:28 -0700 (Wed, 20 Feb 2013)</p><p>BRANCH COMMIT - land ice<br>
<br>
Adding a calving 'law'.<br>
<br>
<br>
This commit just includes a single calving law that calves any floating ice less than a critical thickness. However, the framework is laid to add other calving laws (though some calving implementation may need to occur at other points in the code).<br>
<br>
The changes made are primarily to 1) the implementation of the calving law and 2) modifications to the CFBC sub-gride marine advance parameterization added in recent commits. Details below.<br>
<br>
(Note: There is a difference in the thickness field when running on 1 and 4 processors of ~1e-11. I determined this was present prior to this commit, so I will deal with it in a followup commit.)<br>
<br>
1. Calving Law<br>
This is all done in land_ice_apply_calving(). Right now the critical thickness calving law is the only option. New namelist parameters have been introduced:<br>
namelist character land_ice_model config_calving_law none<br>
% valid options are none, critical-thickness, ocean-kill<br>
namelist real land_ice_model config_calving_critical_thickness 200.0<br>
<br>
2. CFBC-related modifications<br>
Most changes to the CFBC code (adjust_marine_boundary_fluxes) are to properly evolve iceArea.<br>
* Adding a new state variable called iceArea that tracks what fraction of a grid cell is covered by ice. It is 0 for non-ice cells and it is equal to areaCell for grounded ice cells and most floating ice cells. It only has a different value on partially-filled marine boundary cells when config_marine_advance==CFBC. In the code, thickness and layerThickness assume that a grid cell's area is entirely covered by ice (cell fills from bottom up). This variable tracks the amount of area covered by ice in marine-advancing cells, which fill from the side. <br>
* Initializing iceArea in mpas_init_block. Initially there are no partially-filled marine advance cells.<br>
* Scaling SMB and BMB on partially filled marine margin cells to account for iceArea.<br>
* Adding iceArea and layerThickness to the CFBC subroutine (adjust_marine_boundary_fluxes)<br>
* Made numUpwindNeighbors and Href local arrays rather than scalars - we now need their values later in the CFBC subroutine.<br>
* Only consider a neighbor to be upwind if there is a net flux from that cell. Before I just checked the surface velocity direction, but with HO/Stokes solvers there is the possibility that not all layers are moving in the same direction.<br>
* Update iceArea, not allowing a cell to be overfilled based on its reference height. This is a large chunk of somewhat confusing code at the end of the subroutine. It has to find cells that were previously empty but now have had ice spill into them.<br>
<br>
3. Other changes:<br>
* Using thickness variable for the sum of layerThicknesses in the vertical remapping scheme - this eliminate the need to initialize layerThickness explicitly in mpas_init_block<br>
* Adding a check in land_ice_vel_solve to make sure the velocity solver solution is consistent with the mask we gave it. If the velocity solver gives nonzero velocities on thin ice edges, there should be a fatal error, but for now just a warning is printed and the velocity on those edges is set to 0. This is a hack until LifeV can be updated.<br>
* Adding normalVelocity and the 3 masks as output fields. These are useful for development/debugging but are not necessary for production runs.<br>
* Added halo update after velocity solve at initial time.<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-19 17:55:19 UTC (rev 2483)
+++ branches/land_ice_projects/implement_core/src/core_land_ice/Registry 2013-02-20 16:10:28 UTC (rev 2484)
@@ -34,7 +34,12 @@
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
+namelist character land_ice_model config_calving_law none
+% valid options are none, critical-thickness, ocean-kill
+namelist real land_ice_model config_calving_critical_thickness 200.0
+% only used with critical-thickness calving law - this is the ice thickness at which ice calves.
+
% 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
%namelist integer land_ice_model config_thickness_adv_order 1
@@ -211,7 +216,9 @@
var persistent real bedTopography ( nCells Time ) 2 ir 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 ir normalVelocity state - -
+var persistent real normalVelocity ( nVertLevels nEdges Time ) 2 iro normalVelocity state - -
+var persistent real iceArea ( nCells Time ) 2 ro iceArea state - -
+% iceArea is the area of a grid cell that is covered by ice. It is only used for sub-grid marine advance but is defined for all cells.
%var persistent real sup_layerThickness ( nVertLevels nCells Time ) 2 iro sup_layerThickness state sup_layerThickness bob
% Tendency variables
@@ -236,10 +243,10 @@
% Masks: calculated diagnostically at each time step
-var persistent integer cellMask ( nCells Time ) 2 r cellMask state - -
+var persistent integer cellMask ( nCells Time ) 2 ro cellMask state - -
% cellMask only needs to be a restart field if config_allow_additional_advance = false (to keep the mask of initial ice extent)
-var persistent integer edgeMask ( nEdges Time ) 2 - edgeMask state - -
-var persistent integer vertexMask ( nVertices Time ) 2 - vertexMask state - -
+var persistent integer edgeMask ( nEdges Time ) 2 o edgeMask state - -
+var persistent integer vertexMask ( nVertices Time ) 2 o vertexMask state - -
% Needed for LifeV solver to know if the region to solve on the block's domain has changed (treated as a logical):
% THIS SHOULD BE CHANGED TO INTEGER TYPE ONCE THE TRUNK IS MERGED TO THIS BRANCH AND field0dinteger IS DEFINED IN framework/mpas_grid_types.F
var persistent real anyVertexMaskChanged ( ) 2 - anyVertexMaskChanged state - -
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-19 17:55:19 UTC (rev 2483)
+++ branches/land_ice_projects/implement_core/src/core_land_ice/mpas_land_ice_mpas_core.F 2013-02-20 16:10:28 UTC (rev 2484)
@@ -1,3 +1,6 @@
+
+#include "mpas_land_ice_mask.inc"
+
module mpas_core
use mpas_framework
@@ -67,11 +70,11 @@
block => domain % blocklist
do while (associated(block))
! Initialize blocks
- call mpas_init_block(block, block % mesh, dt)
+ call mpas_init_block(block, block % mesh, dt, domain % dminfo)
! update halos on velocity
- call mpas_dmpar_exch_halo_field2d_real(domain % dminfo, block % state % time_levs(1) % state % normalVelocity % array, &
- block % mesh % nVertLevels, block % mesh % nEdges, &
- block % parinfo % edgesToSend, block % parinfo % edgesToRecv)
+ !call mpas_dmpar_exch_halo_field2d_real(domain % dminfo, block % state % time_levs(1) % state % normalVelocity % array, &
+ ! block % mesh % nVertLevels, block % mesh % nEdges, &
+ ! block % parinfo % edgesToSend, block % parinfo % edgesToRecv)
! Assign initial time stamp
block % state % time_levs(1) % state % xtime % scalar = startTimeStamp
@@ -158,7 +161,7 @@
end subroutine simulation_clock_init
- subroutine mpas_init_block(block, mesh, dt)
+ subroutine mpas_init_block(block, mesh, dt, dminfo)
use mpas_grid_types
use mpas_rbf_interpolation
@@ -175,9 +178,9 @@
type (block_type), intent(inout) :: block
type (mesh_type), intent(inout) :: mesh
real (kind=RKIND), intent(in) :: dt
+ type (dm_info), intent(in) :: dminfo
+
type (MPAS_Time_Type) :: currTime
-
- type (dm_info) :: dminfo
integer :: err, err_tmp
type (state_type), pointer :: state
@@ -219,17 +222,16 @@
call mpas_timer_start("initial state calculation")
- ! Need to define initial layerThickness because the vertical remapping that occurs in
- ! land_ice_diagnostic_solve() now uses the old layerThickness values (intead of thickness)
- do iCell = 1, mesh % nCells
- do iLevel = 1, mesh % nVertLevels
- state % layerThickness % array(iLevel,iCell) = mesh % layerThicknessFractions % array(iLevel) &
- * state % thickness % array(iCell)
- end do
- end do
-
call land_ice_calculate_mask_init(mesh, state, err)
call land_ice_diagnostic_solve(mesh, state, err)
+
+ ! Initialize iceArea
+ where (MASK_IS_ICE(state % cellMask % array))
+ state % iceArea % array = mesh % areaCell % array
+ elsewhere
+ state % iceArea % array = 0.0_RKIND
+ end where
+
! Set the initial flag to 1 so LifeV will calculate its mesh information the first time
state % anyVertexMaskChanged % scalar = 1
call mpas_timer_stop("initial state calculation")
@@ -239,12 +241,17 @@
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)
+ ! TODO
+ ! update halos on velocity
+ call mpas_dmpar_exch_halo_field2d_real(dminfo, state % normalVelocity % array, &
+ mesh % nVertLevels, mesh % nEdges, &
+ block % parinfo % edgesToSend, block % parinfo % edgesToRecv)
+ !call mpas_dmpar_exch_halo_field2d_real(domain % dminfo, block % state % time_levs(1) % state % normalVelocity % array, &
+ ! block % mesh % nVertLevels, block % 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, &
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-19 17:55:19 UTC (rev 2483)
+++ branches/land_ice_projects/implement_core/src/core_land_ice/mpas_land_ice_tendency.F 2013-02-20 16:10:28 UTC (rev 2484)
@@ -42,7 +42,8 @@
public :: land_ice_tendency_h, &
land_ice_tendency_tracers, &
land_ice_diagnostic_solve, &
- land_ice_diagnostic_solve_using_velocity
+ land_ice_diagnostic_solve_using_velocity, &
+ land_ice_apply_calving
!--------------------------------------------------------------------
!
@@ -115,9 +116,17 @@
!
!-----------------------------------------------------------------
integer nVertLevels
+ real (kind=RKIND), dimension(:), pointer :: iceArea, areaCell, marineBasalMassBal, sfcMassBal
+ integer, dimension(:), pointer :: cellMask
nVertLevels = mesh % nVertLevels
+ areaCell => mesh % areaCell % array
+ sfcMassBal => mesh % sfcMassBal % array
+ marineBasalMassBal => mesh % marineBasalMassBal % array
+ iceArea => state % iceArea % array
+ cellMask => state % cellMask % array
+
! 0 tendency
layerThickness_tend = 0.0_RKIND
@@ -133,8 +142,8 @@
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)
+ call adjust_marine_boundary_fluxes(mesh, state % thickness % array, state % layerThickness % array, state % layerThicknessEdge % array, &
+ state % cellMask % array, state % bedTopography % array, dt, layerThickness_tend, state % normalVelocity % array, iceArea, err)
!case ('FCT') !===================================================
! MJH TEMP TRYING IT WITH THICKNESS
@@ -143,6 +152,12 @@
! call mpas_ocn_tracer_advection_tend(stateOld % sup_thickness % array, uh , &
! wTop, stateOld % sup_thickness % array(1,:,:), stateOld % sup_thickness % array(1,:,:), &
! dt/SecondsInYear , mesh, 0.0*layerThickness_tend, block % tend % sup_thickness % array)
+
+ ! Doug suggested actually doing it this way:
+ ! call mpas_ocn_tracer_advection_tend(stateOld % sup_thickness % array, normalVelocity , &
+ ! wTop, stateOld % sup_thickness % array(1,:,:) * 0.0 + 1.0 , stateOld % sup_thickness % array(1,:,:) * 0.0 + 1.0, &
+ ! dt/SecondsInYear , mesh, 0.0*layerThickness_tend, block % tend % sup_thickness % array)
+
! where (stateOld % sup_thickness % array .gt. 1.0e-9)
! block % tend % sup_thickness % array = block % tend % sup_thickness % array / stateOld % sup_thickness % array
! else where
@@ -164,20 +179,39 @@
case ('None') !===================================================
! Do nothing - don't add the MB
case default
+
+
+ ! Make some potential adjustments to SMB and BMB before applying them.
+ ! It's ok to overwrite the values with 0's here, because each time step
+ ! we get a fresh copy of the array from the annual_forcing subroutine.
+ ! 1. make adjustments for where the ice is grounded and floating.
+ ! TODO: more complicated treatment at GL?
+ where ( MASK_IS_GROUNDED(cellMask) )
+ ! Apply marineBasalMassBal to floating ice only.
+ marineBasalMassBal = 0.0_RKIND
+ elsewhere ( MASK_IS_FLOATING(cellMask) )
+ ! Apply groundedBasalMassBal to grounded ice only.
+ ! < PLACEHOLDER >
+ end where
+ ! 2. Make adjustments for partially filled cells if CFBC is on.
+ if (config_marine_advance == 'CFBC') then
+ where ( MASK_IS_FLOATING(cellMask) .and. MASK_IS_MARGIN(cellMask) )
+ ! Adjust the SMB and BMB for partially filled cells to give the proper volume flux
+ sfcMassBal = sfcMassBal * iceArea / areaCell
+ marineBasalMassBal = marineBasalMassBal * iceArea / areaCell
+ else where (MASK_IS_NOT_ICE(cellMask) )
+ ! We don't allow a positive BMB where ice is not already present.
+ mesh % marineBasalMassBal % array = 0.0_RKIND
+ end where
+ end if
+
! Add surface mass balance to tendency
- ! TODO: Need to decide where to put this for layer by layer thickness tendency, and how to deal with negative SMB that eliminates top layer or all ice (check for negative thickness?)
- layerThickness_tend(1,:) = layerThickness_tend(1,:) + mesh % sfcMassBal % array ! (tendency in meters per year)
+ ! TODO: Need to decide how to deal with negative SMB that eliminates top layer or all ice (check for negative thickness?)
+ layerThickness_tend(1,:) = layerThickness_tend(1,:) + sfcMassBal ! (tendency in meters per year)
! TODO THIS MIGHT RESULT IN NEGATIVE LAYER THICKNESS!
! Add basal mass balance to tendency
- ! Apply marineBasalMassBal to floating ice only. Apply groundedBasalMassBal to grounded ice only.
- ! TODO: more complicated treatment at GL?
- ! Assign 0.0 values to the marineBasalMassBal where the ice is grounded.
- ! It's ok to overwrite the values with 0's here, because each time step
- ! we get a fresh copy of the array from the annual_forcing subroutine.
- where ( MASK_IS_GROUNDED(state % cellMask % array) )
- mesh % marineBasalMassBal % array = 0.0_RKIND
- end where
+ ! TODO: Need to decide how to deal with negative BMB that eliminates top layer or all ice (check for negative thickness?)
layerThickness_tend(nVertLevels,:) = layerThickness_tend(nVertLevels,:) &
+ mesh % marineBasalMassBal % array ! (tendency in meters per year)
! TODO Add in grounded ice basal mass balance once temperature diffusion is calculated
@@ -439,7 +473,7 @@
! as always, upper surface is the lower surface plus the thickness
upperSurface(:) = lowerSurface(:) + thickness(:)
-
+ ! Do vertical remapping of thickness and tracers
call vertical_remap(layerThickness, thickness, state % tracers % array, mesh, err)
end subroutine land_ice_diagnostic_solve
@@ -535,6 +569,153 @@
end subroutine land_ice_diagnostic_solve_using_velocity
+
+!***********************************************************************
+!
+! subroutine land_ice_apply_calving
+!
+!> \brief Applies a calving 'law' to any marine-terminating ice
+!> \author Matt Hoffman
+!> \date 19 February 2013
+!> \version SVN:$Id$
+!> \details
+!> This routine is a driver for applying various calving laws to any
+!> marine-terminating ice.
+!
+!-----------------------------------------------------------------------
+ subroutine land_ice_apply_calving(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
+ !
+ !-----------------------------------------------------------------
+ integer, dimension(:), pointer :: cellMask, nEdgesOnCell
+ integer, dimension(:,:), pointer :: cellsOnCell
+ real (kind=RKIND), dimension(:), pointer :: thickness, iceArea, areaCell
+ real (kind=RKIND), dimension(:,:), pointer :: layerThickness
+ real (kind=RKIND), dimension(:,:,:), pointer :: tracers
+ integer :: iCell, nCells, numCalvedCells, nThickEdges, iEdge
+ real (kind=RKIND) :: physicalThickness
+
+ areaCell => mesh % areaCell % array
+ nEdgesOnCell => mesh % nEdgesOnCell % array
+ cellsOnCell => mesh % cellsOnCell % array
+
+ cellMask => state % cellMask % array
+ thickness => state % thickness % array
+ iceArea => state % iceArea % array
+ layerThickness => state % layerThickness % array
+ tracers => state % tracers % array
+
+ nCells = mesh % nCells
+
+ select case (config_calving_law)
+
+ !====================================
+ !case ('ocean-kill')
+
+ !====================================
+ case ('critical-thickness')
+ ! This case loops around the floating ice margin, removing any ice below
+ ! a critical thickness. Multiple loops are needed until the entire region
+ ! of too-thin ice has been found.
+ ! This could be moved to its own subroutine...
+ numCalvedCells = 0
+
+ ! Loop across the entire domain multiple times until the inward progression of too-thin ice has been totally mapped
+ do
+
+ ! Update the mask before we look for the margin, since the margin may have advanced.
+ call land_ice_calculate_mask(mesh, state, err)
+
+ do iCell = 1, nCells
+ if ( (MASK_IS_FLOATING(cellMask(iCell))) .and. (MASK_IS_MARGIN(cellMask(iCell))) ) then
+ ! Convert to physical thickness - only needed if CFBC is on but this is always defined.
+ physicalThickness = thickness(iCell) * areaCell(iCell) / iceArea(iCell) ! iceArea should always be > 0 since we are only checking on ice cells.
+ if (iceArea(iCell) == 0.0) then
+ write(0,*) 'iceArea is 0 on a cell with ice!'
+ err = 1
+ endif
+ if (physicalThickness <= config_calving_critical_thickness) then
+ numCalvedCells = numCalvedCells + 1
+ ! TODO tally a sum of mass lost to calving
+ ! Make this cell ice-free and set tracers there to 0
+ thickness(iCell) = 0.0_RKIND
+ layerThickness(:,iCell) = 0.0_RKIND
+ tracers(:,:,iCell) = 0.0_RKIND
+ iceArea(iCell) = 0.0_RKIND
+ end if
+ end if
+ end do
+
+ if (numCalvedCells == 0) then
+ ! if we haven't calved any cells this time though the loop, then we have mapped the extent of the too-thin ice
+ exit
+ else
+ numCalvedCells = 0
+ end if
+
+ end do ! start the search again
+
+ ! clean up - remove any hanging pieces that will not be dynamically active anymore. The mask should already be updated with current ice extent so we don't need to do it again.
+ do iCell = 1, nCells
+ if ( (MASK_IS_FLOATING(cellMask(iCell))) .and. (MASK_IS_MARGIN(cellMask(iCell))) ) then
+ nThickEdges = 0
+ do iEdge = 1, nEdgesOnCell(iCell)
+ if ( MASK_IS_THICK_ICE(cellMask(cellsOnCell(iEdge, iCell))) ) then
+ nThickEdges = nThickEdges + 1
+ endif
+ end do
+ if (nThickEdges == 0) then
+ ! Make this cell ice-free and set tracers there to 0
+ thickness(iCell) = 0.0_RKIND
+ layerThickness(:,iCell) = 0.0_RKIND
+ tracers(:,:,iCell) = 0.0_RKIND
+ iceArea(iCell) = 0.0_RKIND
+ end if
+ end if
+ end do
+
+ ! Update the mask one last time (probably not needed because land_ice_diagnostic solve is coming up)
+ call land_ice_calculate_mask(mesh, state, err)
+
+ !====================================
+ !case ('eigencalving')
+
+ end select
+
+ end subroutine land_ice_apply_calving
+
+
+
+
! private subroutines================================================
@@ -763,9 +944,13 @@
CellUpwind = cell2
CellDownwind = cell1
endif
- maxAllowableDt = (0.5_RKIND * dcEdge(iEdge)) / (abs(normalVelocity(k, iEdge)) + 1.0e-36_RKIND) ! in years
- if ( (dt/SecondsInYear) .gt. maxAllowableDt ) then
- !write(0,*) 'CFL violation at level, edge', k, iEdge
+ if (abs(normalVelocity(k, iEdge)) > 0.0_RKIND) then
+ maxAllowableDt = (0.5_RKIND * dcEdge(iEdge)) / abs(normalVelocity(k, iEdge)) ! in years
+ else
+ maxAllowableDt = 1.0e36_RKIND
+ endif
+ if ( maxAllowableDt < (dt/SecondsInYear) ) then
+ write(0,*) 'CFL violation at level, edge', k, iEdge
err = err + 1
endif
MinOfMaxAllowableDt = min(MinOfMaxAllowableDt, maxAllowableDt)
@@ -778,6 +963,7 @@
if (abs(VelSignSum) .ne. nVertLevels) then
write (6,*) 'Unusual occurrence: Not all layers are moving the same direction on edge ', iEdge
endif
+
end do
if (err .gt. 0) then
@@ -832,7 +1018,7 @@
!
!-----------------------------------------------------------------------
- subroutine adjust_marine_boundary_fluxes(mesh, thickness, layerThicknessEdge, cellMask, bedTopography, deltat, layerThickness_tend, normalVelocity, err)!{{{
+ subroutine adjust_marine_boundary_fluxes(mesh, thickness, layerThickness, layerThicknessEdge, cellMask, bedTopography, deltat, layerThickness_tend, normalVelocity, iceArea, err)!{{{
!-----------------------------------------------------------------
!
@@ -847,6 +1033,9 @@
thickness !< Input: thickness of each cell
real (kind=RKIND), dimension(:,:), intent(in) :: &
+ layerThickness !< Input: thickness of each layer on cells
+
+ real (kind=RKIND), dimension(:,:), intent(in) :: &
layerThicknessEdge !< Input: thickness of each layer on edges
integer, dimension(:), intent(in) :: &
@@ -864,12 +1053,20 @@
!-----------------------------------------------------------------
real (kind=RKIND), dimension(:,:), intent(inout) :: &
- layerThickness_tend !< Input: layer thickness tendency; modified on advancing marine margins
+ layerThickness_tend !< Input/Output: layer thickness tendency; modified on advancing marine margins
real (kind=RKIND), dimension(:,:), intent(inout) :: &
- normalVelocity !< Input: velocity of edges; modified on advancing marine margins
+ normalVelocity !< Input/Output: 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.
+ real (kind=RKIND), dimension(:), intent(inout) :: &
+ iceArea !< Input/Output: area of grid cell covered by ice
+ !< Note: Right now this subroutine modified iceArea directly rather than
+ !< calculating a tendency to be applied later. Therefore this iceArea
+ !< time level needs to be moved to the other time level when time stepping
+ !< actually occurs. This is a kind of confusing implementation and probably
+ !< should be changed.
+
!-----------------------------------------------------------------
!
! output variables
@@ -889,9 +1086,10 @@
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
+ integer :: iCell, nCells, iEdge, neighborIndex, edgeIndex, k, nVertLevels, jEdge, kEdge, jCell
! stuff for making calculations
- real (kind=RKIND) :: numUpwindNeighbors, sumUpwindThickness, numOceanNeighbors, Href, projectedThickness, ExtraVolume, tendencyLayerVolumeRateAdjustment, VelSign, minHref, maxHref
+ real (kind=RKIND) :: sumUpwindThickness, numOceanNeighbors, projectedThickness, ExtraVolume, tendencyLayerVolumeRateAdjustment, VelSign
+ real (kind=RKIND), dimension(:), allocatable :: numUpwindNeighbors, Href ! We need these two to be local arrays
! Only execute the contents of this subroutine if CFBC is turned on.
@@ -900,7 +1098,6 @@
!write(0,*), 'Doing CFBC'
! === Initialize variables and pointers ===
-
nCells = mesh % nCells
nVertLevels = mesh % nVertLevels
@@ -911,9 +1108,10 @@
edgesOnCell => mesh % edgesOnCell % array
cellsOnEdge => mesh % cellsOnEdge % array
- minHref = 99999.0
- maxHref = -1.0
+ allocate(numUpwindNeighbors(nCells+1))
+ allocate(Href(nCells+1))
+
! === 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.
@@ -922,20 +1120,18 @@
! 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
+ ! === 1. Calculate the reference height from the edges feeding this cell ===
+ numUpwindNeighbors(iCell) = 0.0_RKIND
+ sumUpwindThickness = 0.0_RKIND
+ numOceanNeighbors = 0.0_RKIND
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
+
+ ! check which edges have net flux into this cell ! + velocity goes from index 1 to 2 in the cellsOnEdge array.
+ ! Note: We check the entire column's flux because we can't assume that all layers are moving in the same direction
+ if ( ((sum(normalVelocity(:,edgeIndex)*layerThicknessEdge(:, edgeIndex)) > 0.0_RKIND) .and. (cellsOnEdge(2, edgeIndex) == iCell)) &
+ .or. ((sum(normalVelocity(:,edgeIndex)*layerThicknessEdge(:, edgeIndex)) < 0.0_RKIND) .and. (cellsOnEdge(1, edgeIndex) == iCell)) ) then
+ numUpwindNeighbors(iCell) = numUpwindNeighbors(iCell) + 1.0_RKIND
sumUpwindThickness = sumUpwindThickness + sum(layerThicknessEdge(:, edgeIndex))
end if
@@ -943,40 +1139,54 @@
! 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
+ numOceanNeighbors = numOceanNeighbors + 1.0_RKIND
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
+ end do ! loop over edges on cell
+
+
+ ! Calculate reference height for this cell.
+ if (numUpwindNeighbors(iCell) == 0.0_RKIND) then
+ ! It's ok if a cell has no thick ice neighbors now. We will wait until we see if any of its partially filled neighbors spill into it before deciding what to do about it. But set the Href to be defined so we don't enter the if-statement below checking for overfilled cells.
+ Href(iCell) = 1.0e36_RKIND
else
- !write(0,*) 'ok: margin cell has x neighbors', numUpwindNeighbors
+ Href(iCell) = sumUpwindThickness / numUpwindNeighbors(iCell)
+ !write(0,*) 'ok: margin cell has x upwind neighbors', numUpwindNeighbors(iCell)
endif
- Href = sumUpwindThickness / numUpwindNeighbors
- !write(0,*) 'Href=', Href
+
+ !write(0,*) 'Href=', Href(iCell)
! 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 ===
+ !if (sum(layerThickness_tend(:,iCell)) < 0.0) then
+ ! write(0,*) 'thickness tend < 0 for a marine margin cell.'
+ ! write(0,*) layerThickness_tend(:,iCell)
+ !endif
+
+
+ ! === 2. 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
+ ! Update iceArea for this cell, not allowing it to overfill
+ iceArea(iCell) = min(projectedThickness, Href(iCell)) * areaCell(iCell) / Href(iCell)
+
+ if (projectedThickness > Href(iCell)) then
+
! === Calculate residual thickness and move forward into open ocean cells ===
- ExtraVolume = (projectedThickness - Href) * areaCell(iCell) ! units = m^3
+ ExtraVolume = (projectedThickness - Href(iCell)) * areaCell(iCell) ! units = m^3
+ ! Make sure we have somewhere to move the ice to.
+ if (numOceanNeighbors<1) then
+ ! TODO Deal with this situation without a fatal error. This situation is unlikely but possible.
+ 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
+
+ ! Find the open ocean neighbors
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
+ ! TODO this could be weighted by the layer thickness for the case of non-uniform sigma layers
tendencyLayerVolumeRateAdjustment = ExtraVolume / numOceanNeighbors / nVertLevels / (deltat / SecondsInYear)
! Add the tendency to the neighbor. (Need to add it to the existing value
@@ -986,7 +1196,7 @@
layerThickness_tend(:,iCell) = layerThickness_tend(:,iCell) - &
tendencyLayerVolumeRateAdjustment / areaCell(iCell)
- ! === Now calculate a consistent velocity on the edges passing mass to the new cells ===
+ ! === Now calculate a consistent velocity on the edges passing mass to the new cells (needed for tracer advection) ===
! 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)
@@ -1007,8 +1217,8 @@
* VelSign
end do
end if
- end do
- end do
+ end do ! kEdge
+ end do ! jEdge
end if ! open ocean neighbor
@@ -1017,10 +1227,79 @@
end if ! floating margin cell
end do ! cell loop
- !write(0,*) 'CFBC: min, max href=', minHref, maxHref ! Useful for debugging.
+ ! === 3. Now that we have completed the loop through cells, we will have new fluxes into previously open-ocean cells OR cells that would have otherwise been dynamically stranded. We need to update iceArea for them.
+ ! We could not do that before as we went, because we need to know how many upwind neighbors they have, information that is complete now.
+ do iCell = 1, nCells ! TODO - this could probably just be on nCellsSolve
+ ! Identify cells the mask still calls open ocean or are stranded marine thin ice cells but will have ice once time is advanced (have a positive thickness tendency)
+ if ( ( ((MASK_IS_NOT_ICE(cellMask(iCell))) .and. (bedTopography(iCell) < config_sealevel)) .or. &
+ (MASK_IS_THIN_ICE(cellMask(iCell)) .and. MASK_IS_FLOATING(cellMask(iCell)) .and. numUpwindNeighbors(iCell) == 0.0_RKIND) ) .and. &
+ (sum(layerThickness_tend(:,iCell))>0.0_RKIND) ) then
+ ! === Calculate the reference height from the edges feeding this cell ===
+ numUpwindNeighbors(iCell) = 0.0_RKIND
+ sumUpwindThickness = 0.0_RKIND
+ 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 that layerThicknessEdge is no longer valid for these edges because they previously had no flux
+
+ ! these edges should have 0 velocity except for edges adjacent to newly filled cells which will always have all layers fluxing into this cell (what we set in the previous step above), so we don't need to check the entire flux, we can just check a single layer's velocity.
+ if ( ( (normalVelocity(1,edgeIndex) > 0.0_RKIND) .and. (cellsOnEdge(2, edgeIndex) == iCell)) &
+ .or. ( (normalVelocity(1,edgeIndex) < 0.0_RKIND) .and. (cellsOnEdge(1, edgeIndex) == iCell)) ) then
+
+
+! 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
+
+ ! Find the neighbor cell index on this edge
+ do jCell = 1, 2
+ if (cellsOnEdge(jCell, edgeIndex) /= iCell) then
+ neighborIndex = cellsOnEdge(jCell, edgeIndex)
+ endif
+ end do
+
+ numUpwindNeighbors(iCell) = numUpwindNeighbors(iCell) + 1.0_RKIND
+ sumUpwindThickness = sumUpwindThickness + Href(neighborIndex) ! Use the neighbor's Href since it was a partially full cell
+ end if
+ end do
+ if (numUpwindNeighbors(iCell) == 0.0_RKIND) then
+ write(0,*) 'ERROR: There are 0 upwind neighbors but this cell has a tendency'
+ err = 1
+ endif
+
+ Href(iCell) = sumUpwindThickness / numUpwindNeighbors(iCell)
+ !write(0,*) 'Href=', Href
+ ! TODO make this a flux averaged reference height
+
+ ! Now that we have the Href for the newly advancing cell, calculate what its new iceArea will be.
+ projectedThickness = sum(layerThickness_tend(:,iCell)) ! There is no other mass besides the tendency
+ ! TODO - do we want this cell to have SMB/BMB applied??? If so it will probably be over a sliver, but is probably worth accounting for anyway.
+ if (projectedThickness > Href(iCell)) then
+ write(0,*) "Error: in CFBC a newly advanced ice shelf margin cell has overfilled from its neighbors' residuals!"
+ endif
+ iceArea(iCell) = min(projectedThickness, Href(iCell)) * areaCell(iCell) / Href(iCell)
+ endif
+ end do ! cell loop
+
+
+ ! One last step. It is possible for a cell to become dynamically stranded during retreat so that it still has some ice in it but velocity is 0 on all edges feeding it. In the CFBC subgrid parameterization, this cell is effectively in front of the calving front, and should therefore be calved.
+ do iCell = 1, nCells
+ if ( numUpwindNeighbors(iCell) == 0.0_RKIND .and. MASK_IS_THIN_ICE(cellMask(iCell)) .and. MASK_IS_FLOATING(cellMask(iCell)) ) then
+ ! These cells have been dynamically stranded beyond the calving front on retreat, and they don't have any overfilled neighbors to save them - therefore calve them off.
+ layerThickness_tend(:,iCell) = -1.0_RKIND * layerThickness(:,iCell) * (deltat / SecondsInYear)
+ iceArea(iCell) = 0.0_RKIND
+ ! Tracers should be eliminated when the ice is eliminated.
+
+ endif
+ end do
+
+
+ deallocate(numUpwindNeighbors)
+ deallocate(Href)
+
end if ! to do CFBC or not
+
end subroutine adjust_marine_boundary_fluxes
@@ -1233,7 +1512,7 @@
! *** Compute new layer thicknesses (layerInterfaceSigma coordinates)
do iCell = 1, nCells
- thisThk = sum(layerThickness(:,iCell))
+ thisThk = thickness(iCell)
do k = 1, nVertLevels
layerThickness(k,iCell) = layerThicknessFractions(k) * thisThk
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 2013-02-19 17:55:19 UTC (rev 2483)
+++ branches/land_ice_projects/implement_core/src/core_land_ice/mpas_land_ice_time_integration_forwardeuler.F 2013-02-20 16:10:28 UTC (rev 2484)
@@ -23,6 +23,7 @@
use mpas_vector_reconstruction
use land_ice_vel, only: land_ice_vel_solve
use land_ice_tendency
+ use land_ice_mask
contains
@@ -115,7 +116,15 @@
call mpas_dmpar_exch_halo_field2d_real(dminfo, layerthickness_tend, &
mesh % nVertLevels, mesh % nCells, &
block % parinfo % cellsToSend, block % parinfo % cellsToRecv)
+
+ ! Also update halos on iceArea since that might be wrong outside of 0-halos. (don't bother unless CFBC is on)
+ if (config_marine_advance == 'CFBC') then
+ call mpas_dmpar_exch_halo_field1d_real(dminfo, stateOld%iceArea%array, &
+ mesh % nCells, &
+ block % parinfo % cellsToSend, block % parinfo % cellsToRecv)
+ endif
+
!print *,' Max abs thickness tend:', maxval(abs(thickness_tend(1:mesh % nCellsSolve)))
block => block % next
@@ -196,6 +205,7 @@
where (thicknessNew < 0.0_RKIND)
mask = 1
thicknessNew = 0.0_RKIND
+ stateNew % iceArea % array = 0.0_RKIND
end where
if (sum(mask) > 0) then
print *,' Cells with negative thickness (set to 0):',sum(mask)
@@ -209,6 +219,22 @@
deallocate(mask)
+ ! Update ice area - it should be 0 for non-ice cells, some fraction for marine margin cells, and areaCell for all other ice cells. The value for marine margin cells is handled above
+ stateNew%iceArea%array = stateOld%iceArea%array
+ call land_ice_calculate_mask(mesh, stateNew, err)
+ where ( MASK_IS_NOT_ICE(stateNew%cellMask%array) )
+ stateNew%iceArea%array = 0.0_RKIND
+ else where ( MASK_IS_ICE(stateNew%cellMask%array) .and. .not.( MASK_IS_MARGIN(stateNew%cellMask%array) .and. MASK_IS_FLOATING(stateNew%cellMask%array) ) )
+ stateNew%iceArea%array = mesh % areaCell % array
+ end where
+
+ ! ice area as c
+! call mpas_dmpar_exch_halo_field1d_real(dminfo, stateNew%iceArea%array, &
+! mesh % nCells, &
+! block % parinfo % cellsToSend, block % parinfo % cellsToRecv)
+
+
+
! Calculate new tracer values =================
if (config_tracer_advection .ne. 'None') then
do iTracer = 1, size(tracersNew, 1)
@@ -222,6 +248,9 @@
end do
endif
+ ! Apply calving after we have updated the new state - TODO Is this the right place?
+ call land_ice_apply_calving(mesh, stateNew, err)
+
call mpas_timer_stop("calc. new prognostic vars")
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 2013-02-19 17:55:19 UTC (rev 2483)
+++ branches/land_ice_projects/implement_core/src/core_land_ice/mpas_land_ice_vel.F 2013-02-20 16:10:28 UTC (rev 2484)
@@ -12,6 +12,8 @@
!
!-----------------------------------------------------------------------
+#include "mpas_land_ice_mask.inc"
+
module land_ice_vel
use mpas_grid_types
@@ -287,9 +289,16 @@
! local variables
!
!-----------------------------------------------------------------
+ integer :: iEdge, nEdges
+ real (kind=RKIND), dimension(:,:), pointer :: normalVelocity
+ integer, dimension(:), pointer :: edgeMask
err = 0
+ nEdges = mesh % nEdges
+ normalVelocity => state % normalVelocity % array
+ edgeMask => state % edgeMask % array
+
select case (config_dycore)
case ('SIA')
call land_ice_sia_solve(mesh, state, err)
@@ -306,6 +315,17 @@
end select
+ do iEdge = 1, nEdges
+ if ( MASK_IS_THIN_ICE(edgeMask(iEdge)) .and. (maxval(abs(normalVelocity(:,iEdge))) /= 0.0_RKIND) ) then
+ err = 1
+ normalVelocity(:,iEdge) = 0.0_RKIND ! this is a hack because the rest of the code requires this, but this condition should really cause a fatal error.
+ endif
+ enddo
+ if (err == 1) then
+ write(0,*) 'Velocity has been calculated on non-dynamic edges. There is a problem with the velocity solver. Velocity on those edges have been set to 0, but this should be a fatal error.'
+ err = 0 ! a hack to let the code continue until this can be fixed in the velocity solver
+ end if
+
!--------------------------------------------------------------------
end subroutine land_ice_vel_solve
</font>
</pre>