<p><b>mhoffman@lanl.gov</b> 2013-02-20 14:16:27 -0700 (Wed, 20 Feb 2013)</p><p>BRANCH COMMIT - land ice<br>
<br>
Modifications to allow exact match between runs with varying numbers of processors. The main thing is restructuring the land_ice_tend_h_layers_fo_upwind() subroutine to loop over cells instead of edges. To facilitate this I have added a edgeSignOnCell variable to mesh, as is done in the ocean core. I have also rearranged some halo updates. I have tried to put them outside of block loops, which is where they should be once we merge the trunk. Note: I am still seeing difference on varying numbers of processors when running LifeV, but not when I run with SIA.<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-20 20:40:41 UTC (rev 2494)
+++ branches/land_ice_projects/implement_core/src/core_land_ice/Registry 2013-02-20 21:16:27 UTC (rev 2495)
@@ -193,6 +193,12 @@
var persistent real layerCenterSigma ( nVertLevels ) 0 - layerCenterSigma mesh - -
var persistent real layerInterfaceSigma ( nVertLevelsPlus1 ) 0 - layerInterfaceSigma mesh - -
+% Sign fields, for openmp and bit reproducibility without branching statements.
+% Right now, only edgeSignOnCell is used, but I am including these others for possible future use
+var persistent integer edgeSignOnCell ( maxEdges nCells ) 0 - edgeSignOnCell mesh - -
+%var persistent integer edgeSignOnVertex ( maxEdges nVertices ) 0 - edgeSignOnVertex mesh - -
+%var persistent integer kiteIndexOnCell ( maxEdges nCells ) 0 - kiteIndexOnCell mesh - -
+
% Boundary Conditions (defined as mesh variables because they are read in once and held constant for the simulation - at least for now)
%
% Note: beta should be moved to diagnostic fields when a process model comes online
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-20 20:40:41 UTC (rev 2494)
+++ branches/land_ice_projects/implement_core/src/core_land_ice/mpas_land_ice_mask.F 2013-02-20 21:16:27 UTC (rev 2495)
@@ -313,6 +313,8 @@
end do
+ ! vertexMask and edgeMask needs halo updates before they can be used. Halo updates should occur outside of block loops.
+
!--------------------------------------------------------------------
end subroutine land_ice_calculate_mask
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-20 20:40:41 UTC (rev 2494)
+++ branches/land_ice_projects/implement_core/src/core_land_ice/mpas_land_ice_mpas_core.F 2013-02-20 21:16:27 UTC (rev 2495)
@@ -195,6 +195,7 @@
! Call init routines ====
call land_ice_setup_vertical_coords(mesh, err_tmp)
err = ior(err, err_tmp)
+ call land_ice_setup_sign_and_index_fields(mesh)
currTime = mpas_get_clock_time(clock, MPAS_NOW, err_tmp)
err = ior(err, err_tmp)
call land_ice_assign_annual_forcing(currTime, mesh, err_tmp)
@@ -215,6 +216,8 @@
call mpas_timer_start("initialize velocity")
call land_ice_vel_block_init(block, err_tmp)
call mpas_timer_stop("initialize velocity")
+
+ call land_ice_calculate_mask_init(mesh, state, err)
err = ior(err, err_tmp)
@@ -222,9 +225,16 @@
call mpas_timer_start("initial state calculation")
- call land_ice_calculate_mask_init(mesh, state, err)
call land_ice_diagnostic_solve(mesh, state, err)
+ ! vertexMask and edgeMask needs halo updates before they can be used. Halo updates should occur outside of block loops.
+ call mpas_dmpar_exch_halo_field1d_integer(dminfo, state % vertexMask % array, &
+ mesh % nVertices, &
+ block % parinfo % verticesToSend, block % parinfo % verticesToRecv)
+ call mpas_dmpar_exch_halo_field1d_integer(dminfo, state % edgeMask % array, &
+ mesh % nEdges, &
+ block % parinfo % edgesToSend, block % parinfo % edgesToRecv)
+
! Initialize iceArea
where (MASK_IS_ICE(state % cellMask % array))
state % iceArea % array = mesh % areaCell % array
@@ -264,6 +274,10 @@
call land_ice_diagnostic_solve_using_velocity(mesh, state, err) ! Some diagnostic variables require velocity to compute
+ call mpas_dmpar_exch_halo_field2d_real(dminfo, state % layerThicknessEdge % array, &
+ mesh % nVertLevels, mesh % nEdges, &
+ block % parinfo % edgesToSend, block % parinfo % edgesToRecv)
+
if(err == 1) then
print *, "An error has occurred in mpas_init_block. Aborting..."
call mpas_dmpar_abort(dminfo)
@@ -531,7 +545,7 @@
!***********************************************************************
!
-! routine and_ice_setup_vertical_coords
+! routine land_ice_setup_vertical_coords
!
!> \brief Initializes vertical coord system
!> \author Matt Hoffman
@@ -611,4 +625,70 @@
end subroutine land_ice_setup_vertical_coords
+ subroutine land_ice_setup_sign_and_index_fields(mesh)!{{{
+
+ type (mesh_type), intent(inout) :: mesh
+
+ integer, dimension(:), pointer :: nEdgesOnCell
+ integer, dimension(:,:), pointer :: edgesOnCell, cellsOnEdge !, edgesOnVertex, cellsOnVertex, verticesOnCell, verticesOnEdge
+ integer, dimension(:,:), pointer :: edgeSignOnCell !, edgeSignOnVertex, kiteIndexOnCell
+
+ integer :: nCells !, nVertices, vertexDegree
+ integer :: iCell, iEdge, iVertex, i, j, k
+
+ nCells = mesh % nCells
+ !nVertices = mesh % nVertices
+ !vertexDegree = mesh % vertexDegree
+
+ nEdgesOnCell => mesh % nEdgesOnCell % array
+ edgesOnCell => mesh % edgeSOnCell % array
+ !edgesOnVertex => mesh % edgesOnVertex % array
+ !cellsOnVertex => mesh % cellsOnVertex % array
+ cellsOnEdge => mesh % cellsOnEdge % array
+ !verticesOnCell => mesh % verticesOnCell % array
+ !verticesOnEdge => mesh % verticesOnEdge % array
+ edgeSignOnCell => mesh % edgeSignOnCell % array
+ !edgeSignOnVertex => mesh % edgeSignOnVertex % array
+ !kiteIndexOnCell => mesh % kiteIndexOnCell % array
+
+ edgeSignOnCell = 0.0_RKIND
+ !edgeSignOnVertex = 0.0_RKIND
+ !kiteIndexOnCell = 0.0_RKIND
+
+ do iCell = 1, nCells
+ do i = 1, nEdgesOnCell(iCell)
+ iEdge = edgesOnCell(i, iCell)
+ !iVertex = verticesOnCell(i, iCell)
+
+ ! Vector points from cell 1 to cell 2
+ if(iCell == cellsOnEdge(1, iEdge)) then
+ edgeSignOnCell(i, iCell) = -1
+ else
+ edgeSignOnCell(i, iCell) = 1
+ end if
+
+ !do j = 1, vertexDegree
+ ! if(cellsOnVertex(j, iVertex) == iCell) then
+ ! kiteIndexOnCell(i, iCell) = j
+ ! end if
+ !end do
+ end do
+ end do
+
+ !do iVertex = 1, nVertices
+ ! do i = 1, vertexDegree
+ ! iEdge = edgesOnVertex(i, iVertex)
+ !
+ ! ! Vector points from vertex 1 to vertex 2
+ ! if(iVertex == verticesOnEdge(1, iEdge)) then
+ ! edgeSignOnVertex(i, iVertex) = -1
+ ! else
+ ! edgeSignOnVertex(i, iVertex) = 1
+ ! end if
+ ! end do
+ !end do
+
+ end subroutine land_ice_setup_sign_and_index_fields!}}}
+
+
end module mpas_core
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-20 20:40:41 UTC (rev 2494)
+++ branches/land_ice_projects/implement_core/src/core_land_ice/mpas_land_ice_tendency.F 2013-02-20 21:16:27 UTC (rev 2495)
@@ -139,7 +139,7 @@
! 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, state % edgeMask % array, layerThickness_tend, dt, allowableDt, err)
+ state % layerThicknessEdge % 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 % layerThickness % array, state % layerThicknessEdge % array, &
@@ -554,8 +554,8 @@
! Calculate h on edges using first order
VelSign = sign(1.0_RKIND, normalVelocity(k, iEdge))
layerThicknessEdge(k,iEdge) = max(VelSign * layerThickness(k, cell1), VelSign * (-1.0_RKIND) * layerThickness(k, cell2))
- ! + velocity goes from index 1 to 2 in the cellsOnEdge array.
- ! Doug does the calculation as: h_edge = max(VelSign, 0.0) * h1 - min(VelSign, 0.0) * h2
+ ! + velocity goes from index 1 to 2 in the cellsOnEdge array.
+ ! Doug does the calculation as: h_edge = max(VelSign, 0.0) * h1 - min(VelSign, 0.0) * h2
end do
! thickness_edge is not currently in registry and not currenly needed. If it is, uncomment the next line
!h_edge = max(thickness(cell1), thickness(cell2))
@@ -566,6 +566,7 @@
endif
! Note: no halo update needed on layerThicknessEdge because normalVelocity has already been halo updated.
+ ! However, the outmost layerThicknessEdge may be wrong if its upwind cell is off this block.
end subroutine land_ice_diagnostic_solve_using_velocity
@@ -860,7 +861,7 @@
!> results.
!
!-----------------------------------------------------------------------
- subroutine land_ice_tend_h_layers_fo_upwind(grid, normalVelocity, layerThickness, edgeMask, tend, dt, MinOfMaxAllowableDt, err)!{{{
+ subroutine land_ice_tend_h_layers_fo_upwind(grid, normalVelocity, layerThicknessEdge, edgeMask, tend, dt, MinOfMaxAllowableDt, err)!{{{
!-----------------------------------------------------------------
!
@@ -875,7 +876,7 @@
normalVelocity !< Input: velocity
real (kind=RKIND), dimension(:,:), intent(in) :: &
- layerThickness !< Input: thickness of each layer
+ layerThicknessEdge !< Input: thickness of each layer on edges
integer, dimension(:), intent(in) :: &
edgeMask !< Input: mask on edges
@@ -906,67 +907,52 @@
!
!-----------------------------------------------------------------
- integer :: iEdge, nEdges, cell1, cell2, nVertLevels, k, CellDownwind, CellUpwind
-
- integer, dimension(:,:), pointer :: cellsOnEdge
-
- real (kind=RKIND) :: flux, VelSign, maxAllowableDt, VelSignSum
+ integer :: nCells, nVertLevels, iEdge, iCell, i, k
+ integer, dimension(:,:), pointer :: edgesOnCell, edgeSignOnCell
+ integer, dimension(:), pointer :: nEdgesOnCell
+ real (kind=RKIND) :: invAreaCell, flux, maxAllowableDt
real (kind=RKIND), dimension(:), pointer :: dvEdge, areaCell, dcEdge
! Only needed for optional check for mass conservation
- integer :: iCell
- real :: tendVolSum
+ !real (kind=RKIND) :: tendVolSum
err = 0
- MinOfMaxAllowableDt = 1.0e36_RKIND
- nEdges = grid % nEdges
+ nCells = grid % nCells
nVertLevels = grid % nVertLevels
-
- cellsOnEdge => grid % cellsOnEdge % array
+ nEdgesOnCell => grid % nEdgesOnCell % array
+ edgesOnCell => grid % edgesOnCell % array
+ edgeSignOnCell => grid % edgeSignOnCell % array
dvEdge => grid % dvEdge % array
dcEdge => grid % dcEdge % array
areaCell => grid % areaCell % array
- do iEdge=1,nEdges
- cell1 = cellsOnEdge(1,iEdge)
- cell2 = cellsOnEdge(2,iEdge)
+ MinOfMaxAllowableDt = 1.0e36_RKIND
- VelSignSum = 0.0
- do k=1, nVertLevels
- ! Don't assume all layers are flowing the same direction, so determine the upwind direction separately for each layer
- VelSign = sign(1.0, normalVelocity(k, iEdge))
- VelSignSum = VelSignSum + VelSign
- if (VelSign .gt. 0.0) then
- CellUpwind = cell1
- CellDownwind = cell2
- else
- CellUpwind = cell2
- CellDownwind = cell1
- endif
- 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)
+ do iCell = 1, nCells
+ invAreaCell = 1.0 / areaCell(iCell)
+ do i = 1, nEdgesOnCell(iCell)
+ iEdge = edgesOnCell(i, iCell)
+ do k = 1, nVertLevels
- 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
- ! Sanity check on ice motion:
- if (abs(VelSignSum) .ne. nVertLevels) then
- write (6,*) 'Unusual occurrence: Not all layers are moving the same direction on edge ', iEdge
+ 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)
+ flux = normalVelocity(k, iEdge) * dvEdge(iEdge) * layerThicknessEdge(k, iEdge)
+ tend(k, iCell) = tend(k, iCell) + edgeSignOnCell(i, iCell) * flux * invAreaCell
+ end do
+ end do
end do
- if (err .gt. 0) then
+ if (err > 0) then
write(0,*) 'CFL violation on this processor on ', err, ' level-edges! Maximum allowable time step (yrs) for this processor is ', MinOfMaxAllowableDt
err = 1
endif
@@ -987,8 +973,6 @@
-
-
!***********************************************************************
!
! subroutine adjust_marine_boundary_fluxes
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-20 20:40:41 UTC (rev 2494)
+++ branches/land_ice_projects/implement_core/src/core_land_ice/mpas_land_ice_time_integration_forwardeuler.F 2013-02-20 21:16:27 UTC (rev 2495)
@@ -71,65 +71,43 @@
dminfo => domain % dminfo
procVertexMaskChanged = 0
+! === Implicit column physics (vertical temperature diffusion) ===========
+ !call ()
+! === Calculate Tendencies ========================
+ call mpas_timer_start("calculate tendencies")
+
+ ! Thickness tendencies
block => domain % blocklist
do while (associated(block))
-! === Initialize pointers ========================
-
- ! \todo: These pointer allocations need to be removed from the block loop before multiple blocks will work!
-
! Mesh information
mesh => block % mesh
- nVertLevels = mesh % nVertLevels
- layerThicknessFractions => mesh % layerThicknessFractions % array
-
- ! State at time n
stateOld => block % state % time_levs(1) % state
- thicknessOld => stateOld % thickness % array
- layerThicknessOld => stateOld % layerThickness % array
- normalVelocityOld => stateOld % normalVelocity % array
- tracersOld => stateOld % tracers % array
- cellMaskOld => stateOld % cellMask % array
-
- ! State at time n+1 (advanced by dt by Forward Euler)
- stateNew => block % state % time_levs(2) % state
- thicknessNew => stateNew % thickness % array
- layerThicknessNew => stateNew % layerThickness % array
- normalVelocityNew => stateNew % normalVelocity % array
- tracersNew => stateNew % tracers % array
-
- ! Tendencies
layerThickness_tend => block % tend % layerThickness % array
- tracer_tendency => block % tend % tracers % array
-
-! === Implicit column physics (vertical temperature diffusion) ===========
- !call ()
-
-
-! === Calculate Tendencies ========================
- call mpas_timer_start("calculate tendencies")
! Calculate thickness tendency using state at time n =========
call land_ice_tendency_h(mesh, stateOld, layerThickness_tend, dt, dminfo, allowableDt, err)
-
- ! update halos on thickness tend (perhaps move to land_ice_tendency_h after we've updated the halo update calls)
- 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
end do
+ ! Now that we have exited the block loop, do any needed halo updates.
+ ! update halos on thickness tend
+ call mpas_dmpar_exch_halo_field2d_real(dminfo, layerthickness_tend, &
+ mesh % nVertLevels, mesh % nCells, &
+ domain % blocklist % parinfo % cellsToSend, domain % blocklist % 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, &
+ domain % blocklist % parinfo % cellsToSend, domain % blocklist % parinfo % cellsToRecv)
+ ! update halos on velocity - this may have been changed
+ call mpas_dmpar_exch_halo_field2d_real(dminfo, stateOld % normalVelocity % array, &
+ mesh % nVertLevels, mesh % nEdges, &
+ domain % blocklist % parinfo % edgesToSend, domain % blocklist % parinfo % edgesToRecv)
+ endif
+
! Determine CFL limit on all procs
call mpas_dmpar_min_real(dminfo, allowableDt, allowableDtMin)
! Determine which processor has the limiting CFL
@@ -146,36 +124,70 @@
call mpas_dmpar_abort(dminfo)
endif
+ ! Tracer tendencies
block => domain % blocklist
do while (associated(block))
+ ! Mesh information
+ mesh => block % mesh
+ ! State at time n
+ stateOld => block % state % time_levs(1) % state
+ ! Tendencies
+ layerThickness_tend => block % tend % layerThickness % array
+ tracer_tendency => block % tend % tracers % array
+
! Calculate tracer tendencies ==========
! There could be a negative layer thickness with SMB turned on!
call land_ice_tendency_tracers(mesh, stateOld, layerThickness_tend, tracer_tendency, dt, dminfo, err)
- if(err == 1) then
- call mpas_dmpar_global_abort("An error has occurred in land_ice_tendency_tracers. Aborting...")
- endif
+ block => block % next
+ end do
- ! update halos on tracer tend (perhaps move to land_ice_tendency_tracers after we've updated the halo update calls)
- select case (config_tracer_advection)
- case ('None') !===================================================
- ! Do nothing - no need to waste time doing a halo update if not advecting tracers! The tendency will be 0 everywhere
- case default
- call mpas_dmpar_exch_halo_field3d_real(dminfo, tracer_tendency, &
- size(tracer_tendency,dim=1), mesh % nVertLevels, mesh % nCells, &
- block % parinfo % cellsToSend, block % parinfo % cellsToRecv)
- end select
+ if (err == 1) then
+ call mpas_dmpar_global_abort("An error has occurred in land_ice_tendency_tracers. Aborting...")
+ endif
- block => block % next
- end do
+ ! update halos on tracer tend
+ select case (config_tracer_advection)
+ case ('None') !===================================================
+ ! Do nothing - no need to waste time doing a halo update if not advecting tracers! The tendency will be 0 everywhere
+ case default
+ call mpas_dmpar_exch_halo_field3d_real(dminfo, tracer_tendency, &
+ size(tracer_tendency,dim=1), mesh % nVertLevels, mesh % nCells, &
+ domain % blocklist % parinfo % cellsToSend, domain % blocklist % parinfo % cellsToRecv)
+ end select
+
call mpas_timer_stop("calculate tendencies")
-
+
+! === Compute new state for prognostic variables ==================================
+! (once implicit column physics are added (i.e. temp diffusion), these calculations will need to be adjusted to apply to the new values as needed)
+ call mpas_timer_start("calc. new prognostic vars")
block => domain % blocklist
do while (associated(block))
-! === Compute new state for prognostic variables ==================================
-! (once implicit column physics are added (i.e. temp diffusion), these calculations will need to be adjusted to apply to the new values as needed)
- call mpas_timer_start("calc. new prognostic vars")
+ ! Mesh information
+ mesh => block % mesh
+ layerThicknessFractions => mesh % layerThicknessFractions % array
+
+ ! State at time n
+ stateOld => block % state % time_levs(1) % state
+ thicknessOld => stateOld % thickness % array
+ layerThicknessOld => stateOld % layerThickness % array
+ normalVelocityOld => stateOld % normalVelocity % array
+ tracersOld => stateOld % tracers % array
+ cellMaskOld => stateOld % cellMask % array
+
+ ! State at time n+1 (advanced by dt by Forward Euler)
+ stateNew => block % state % time_levs(2) % state
+ thicknessNew => stateNew % thickness % array
+ layerThicknessNew => stateNew % layerThickness % array
+ normalVelocityNew => stateNew % normalVelocity % array
+ tracersNew => stateNew % tracers % array
+
+ ! Tendencies
+ layerThickness_tend => block % tend % layerThickness % array
+ tracer_tendency => block % tend % tracers % array
+
+
! Update thickness ======================
! Commented out usage for advecting thickness as a column
@@ -228,13 +240,7 @@
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)
@@ -250,29 +256,43 @@
! Apply calving after we have updated the new state - TODO Is this the right place?
call land_ice_apply_calving(mesh, stateNew, err)
+ block => block % next
+ end do
+ call mpas_timer_stop("calc. new prognostic vars")
- call mpas_timer_stop("calc. new prognostic vars")
-
! === Calculate diagnostic variables for new state =====================
- call mpas_timer_start("calc. diagnostic vars except vel")
+ call mpas_timer_start("calc. diagnostic vars except vel")
+ block => domain % blocklist
+ do while (associated(block))
+ ! Mesh information
+ mesh => block % mesh
+ ! State at time n+1 (advanced by dt by Forward Euler)
+ stateNew => block % state % time_levs(2) % state
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, &
+ block => block % next
+ end do
+
+ ! 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, &
+ domain % blocklist % parinfo % verticesToSend, domain % blocklist % parinfo % verticesToRecv)
+ call mpas_dmpar_exch_halo_field1d_integer(dminfo, stateNew % edgeMask % array, &
mesh % nEdges, &
- block % parinfo % edgesToSend, block % parinfo % edgesToRecv)
+ domain % blocklist % parinfo % edgesToSend, domain % blocklist % parinfo % edgesToRecv)
+ block => domain % blocklist
+ do while (associated(block))
+ stateNew => block % state % time_levs(2) % state
+ stateOld => block % state % time_levs(1) % state
! 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
+ if ( sum(stateNew % vertexMask % array - stateOld % vertexMask % array) /= 0 ) then
blockVertexMaskChanged = 1
else
blockVertexMaskChanged = 0
@@ -282,18 +302,31 @@
! Determine if any blocks on this processor had a change to the vertex mask
procVertexMaskChanged = max(procVertexMaskChanged, blockVertexMaskChanged)
- call mpas_timer_stop("calc. diagnostic vars except vel")
-
-
block => block % next
end do
! Determine if the vertex mask has changed on any processor (need to exit the block loop to do so)
call mpas_dmpar_max_int(dminfo, procVertexMaskChanged, anyVertexMaskChanged)
+ call mpas_timer_stop("calc. diagnostic vars except vel")
+
+
+
call mpas_timer_start("velocity solve")
+
+ ! TODO Once multiple blocks are supported, this section will need to change.
+ ! LifeV does not support multiple blocks but the MPAS SIA could.
block => domain % blocklist
do while (associated(block))
+ ! Mesh information
+ mesh => block % mesh
+ nVertLevels = mesh % nVertLevels
+ ! State at time n
+ stateOld => block % state % time_levs(1) % state
+ normalVelocityOld => stateOld % normalVelocity % array
+ ! State at time n+1 (advanced by dt by Forward Euler)
+ stateNew => block % state % time_levs(2) % state
+ normalVelocityNew => stateNew % normalVelocity % array
! Assign the vertex-changed flag to each block
stateNew % anyVertexMaskChanged % scalar = anyVertexMaskChanged
@@ -304,11 +337,18 @@
normalVelocityNew = normalVelocityOld
call land_ice_vel_solve(mesh, stateNew, err) ! ****** Calculate Velocity ******
- ! update halos on velocity
- call mpas_dmpar_exch_halo_field2d_real(dminfo, normalVelocityNew, &
- nVertLevels, mesh % nEdges, &
- block % parinfo % edgesToSend, block % parinfo % edgesToRecv)
+ block => block % next
+ end do
+ ! update halos on velocity
+ call mpas_dmpar_exch_halo_field2d_real(dminfo, normalVelocityNew, &
+ mesh % nVertLevels, mesh % nEdges, &
+ domain % blocklist % parinfo % edgesToSend, domain % blocklist % parinfo % edgesToRecv)
+
+ block => domain % blocklist
+ do while (associated(block))
+ mesh => block % mesh
+ stateNew => block % state % time_levs(2) % state
! Calculate reconstructed velocities for this time step - do this after velocity halo update in case velocities on the 1-halo edge are wrong (depends on velocity solver)
call mpas_reconstruct(mesh, stateNew % normalVelocity % array,&
stateNew % uReconstructX % array, &
@@ -326,9 +366,14 @@
call mpas_timer_start("calc. diagnostic vars except vel")
block => domain % blocklist
do while (associated(block))
- call land_ice_diagnostic_solve_using_velocity(mesh, stateNew, err) ! Some diagnostic variables require velocity to compute
+ call land_ice_diagnostic_solve_using_velocity(block % mesh, block % state % time_levs(2) % state, err) ! Some diagnostic variables require velocity to compute
block => block % next
end do
+
+ call mpas_dmpar_exch_halo_field2d_real(dminfo, domain % blocklist % state % time_levs(2) % state % layerThicknessEdge % array, &
+ domain % blocklist % mesh % nVertLevels, domain % blocklist % mesh % nEdges, &
+ domain % blocklist % parinfo % edgesToSend, domain % blocklist % parinfo % edgesToRecv)
+
call mpas_timer_stop("calc. diagnostic vars except vel")
! === Cleanup & Misc. =============================
</font>
</pre>