<p><b>mhoffman@lanl.gov</b> 2013-02-04 15:46:17 -0700 (Mon, 04 Feb 2013)</p><p>BRANCH COMMIT - land ice<br>
Some bug fixes and adjustments for improving mass conservation. Note, for the floating dome test case I have been using, mass is still not conserved past the 13th significant digit.<br>
<br>
1) Check on init that layerThicknessFractions exactly sum to 1.0 and fix if they don't.<br>
2) Use, e.g., 1.0_RKIND throughout the code instead of 1.0. Changing this did not seem to make a difference, but since I have done it, there is no reason to revert it. Doug J. says XLF compilers require this (but that's less of an issue with the shutdown of Bluefire).<br>
3) Check during vertical remapping in land_ice_diagnostic_solve for mass conservation and fix if it doesn't.<br>
<br>
Also changed the layer thickness tendency calculation loop (subroutine land_ice_tend_h_layers_fo_upwind) to occur over all edges, not just THICK_ICE edges. All other edges SHOULD have 0 velocity so calculating tendencies on them is fine. (This should have not affect on mass conservation.)<br>
<br>
Also adding sumCellArea to the diagnostic output file and adjusting the write statement of diagnostics to only print significant digits.<br>
</p><hr noshade><pre><font color="gray">Modified: branches/land_ice_projects/implement_core/src/core_land_ice/mpas_land_ice_global_diagnostics.F
===================================================================
--- branches/land_ice_projects/implement_core/src/core_land_ice/mpas_land_ice_global_diagnostics.F 2013-02-04 21:59:39 UTC (rev 2425)
+++ branches/land_ice_projects/implement_core/src/core_land_ice/mpas_land_ice_global_diagnostics.F 2013-02-04 22:46:17 UTC (rev 2426)
@@ -42,6 +42,8 @@
integer, intent(in) :: timeIndex
real (kind=RKIND), intent(in) :: dt
+ character(LEN=20) :: realPrecisionStr
+
integer :: nVertLevels, nCellsSolve, nEdgesSolve, nVerticesSolve, nCellsGlobal, nEdgesGlobal, nVerticesGlobal, iTracer
integer :: nCells
@@ -139,12 +141,12 @@
!!! allocate(h_s_edge(nEdgesSOlve))
- cellVolume = 0
- refAreaWeightedSurfaceHeight = 0
- refAreaWeightedSurfaceHeight_edge = 0
- vertexVolume = 0
- cellArea = 0
- averageThickness = 0
+ cellVolume = 0.0_RKIND
+ refAreaWeightedSurfaceHeight = 0.0_RKIND
+ refAreaWeightedSurfaceHeight_edge = 0.0_RKIND
+ vertexVolume = 0.0_RKIND
+ cellArea = 0.0_RKIND
+ averageThickness = 0.0_RKIND
!volumeWeightedPotentialVorticity = 0
!volumeWeightedPotentialEnstrophy = 0
!volumeWeightedKineticEnergy = 0
@@ -283,7 +285,10 @@
! write(fileID,'(1i0, 100es24.16)') timeIndex, timeIndex*dt, globalFluidThickness, globalPotentialVorticity, globalPotentialEnstrophy, &
! globalEnergy, globalCoriolisEnergyTendency, globalKineticEnergyTendency+globalPotentialEnergyTendency, &
! globalKineticEnergy, globalPotentialEnergy
- write(fileID,'(1i0, 100es24.16)') timeIndex, timeIndex*dt, sumCellVolume
+
+ ! Get the precision being used so we don't write more digits than are meaningful
+ write(realPrecisionStr,*) PRECISION(sumCellVolume)-1
+ write(fileID,'(1i0, 100es24.'// ADJUSTL(realPrecisionStr) // ')') timeIndex, timeIndex*dt, sumCellArea, sumCellVolume
close(fileID)
end if
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-04 21:59:39 UTC (rev 2425)
+++ branches/land_ice_projects/implement_core/src/core_land_ice/mpas_land_ice_mpas_core.F 2013-02-04 22:46:17 UTC (rev 2426)
@@ -40,7 +40,7 @@
! Initialize core
!
! If dt in years is supplied, use it. Otherwise use dt in seconds
- if (config_dt_years .eq. 0.0) then
+ if (config_dt_years .eq. 0.0_RKIND) then
dt = config_dt_seconds
else
dt = config_dt_years * SecondsInYear
@@ -172,12 +172,16 @@
type (state_type), pointer :: state
+ integer :: iCell, iLevel
+
+
err = 0
state => block % state % time_levs(1) % state ! initial state
! Call init routines ====
- call land_ice_setup_vertical_coords(mesh)
+ call land_ice_setup_vertical_coords(mesh, err_tmp)
+ err = ior(err, err_tmp)
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)
@@ -202,7 +206,18 @@
! Initialize state ====
+
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)
! Set the initial flag to 1 so LifeV will calculate its mesh information the first time
@@ -512,7 +527,7 @@
!
!-----------------------------------------------------------------------
- subroutine land_ice_setup_vertical_coords(mesh)
+ subroutine land_ice_setup_vertical_coords(mesh, err)
!-----------------------------------------------------------------
!
@@ -534,6 +549,7 @@
!
!-----------------------------------------------------------------
+ integer, intent(out) :: err
!-----------------------------------------------------------------
!
@@ -542,26 +558,41 @@
!-----------------------------------------------------------------
integer :: nVertLevels, k
+ real (kind=RKIND) :: fractionTotal
real (kind=RKIND), dimension(:), pointer :: layerThicknessFractions, layerCenterSigma, layerInterfaceSigma
nVertLevels = mesh % nVertLevels
! layerThicknessFractions is provided by input
layerThicknessFractions => mesh % layerThicknessFractions % array
+
+ ! Check that layerThicknessFractions are valid
+ ! TODO - switch to having the user input the sigma levels instead???
+ fractionTotal = sum(layerThicknessFractions)
+ if (fractionTotal /= 1.0_RKIND) then
+ if (abs(fractionTotal - 1.0_RKIND) > 0.001_RKIND) then
+ write(0,*) 'The sum of layerThicknessFractions is different from 1.0 by more than 0.001. Aborting.'
+ err = 1
+ end if
+ write (6,*), 'Adjusting upper layerThicknessFrac by small amount because sum of layerThicknessFractions is slightly different from 1.0.'
+ ! TODO - distribute the residual amongst all layers (and then put the residual of that in a single layer
+ layerThicknessFractions(1) = layerThicknessFractions(1) - (fractionTotal - 1.0_RKIND)
+ endif
+
+
! layerCenterSigma is the fractional vertical position (0-1) of each layer center, with 0.0 at the ice surface and 1.0 at the ice bed
layerCenterSigma => mesh % layerCenterSigma % array
! layerInterfaceSigma is the fractional vertical position (0-1) of each layer interface, with 0.0 at the ice surface and 1.0 at the ice bed. Interface 1 is the surface, interface 2 is between layers 1 and 2, etc., and interface nVertLevels+1 is the bed.
layerInterfaceSigma => mesh % layerInterfaceSigma % array
- layerCenterSigma(1) = 0.5 * layerThicknessFractions(1)
- layerInterfaceSigma(1) = 0.0
+ layerCenterSigma(1) = 0.5_RKIND * layerThicknessFractions(1)
+ layerInterfaceSigma(1) = 0.0_RKIND
do k = 2, nVertLevels
- layerCenterSigma(k) = layerCenterSigma(k-1) + 0.5 * layerThicknessFractions(k-1) &
- + 0.5 * layerThicknessFractions(k)
+ layerCenterSigma(k) = layerCenterSigma(k-1) + 0.5_RKIND * layerThicknessFractions(k-1) &
+ + 0.5_RKIND * layerThicknessFractions(k)
layerInterfaceSigma(k) = layerInterfaceSigma(k-1) + layerThicknessFractions(k-1)
end do
- layerInterfaceSigma(nVertLevels+1) = 1.0
+ layerInterfaceSigma(nVertLevels+1) = 1.0_RKIND
-
end subroutine land_ice_setup_vertical_coords
Modified: branches/land_ice_projects/implement_core/src/core_land_ice/mpas_land_ice_sia.F
===================================================================
--- branches/land_ice_projects/implement_core/src/core_land_ice/mpas_land_ice_sia.F 2013-02-04 21:59:39 UTC (rev 2425)
+++ branches/land_ice_projects/implement_core/src/core_land_ice/mpas_land_ice_sia.F 2013-02-04 22:46:17 UTC (rev 2426)
@@ -289,10 +289,10 @@
rhoi = config_ice_density
- basalVelocity = 0.0 ! Assume no sliding
+ basalVelocity = 0.0_RKIND ! Assume no sliding
! Calculate ratefactor (A) at edge - This should be calculated external to this subroutine and as a function of temperature
- ratefactor = 1.0e-16
+ ratefactor = 1.0e-16_RKIND
! Loop over edges
do iEdge = 1, nEdges
@@ -305,14 +305,14 @@
! This could/should be calculated externally to this subroutine
slopeOnEdge = (thickness(cell1) - thickness(cell2) ) / dcEdge(iEdge)
! Calculate thickness on edge - 2nd order
- thicknessEdge = (thickness(cell1) + thickness(cell2) ) * 0.5
+ thicknessEdge = (thickness(cell1) + thickness(cell2) ) * 0.5_RKIND
! Loop over layers
do iLevel = 1, nVertLevels
! Determine the height of each layer above the bed
- layerCenterHeightOnEdge = thicknessEdge * (1.0 - layerCenterSigma(iLevel) )
+ layerCenterHeightOnEdge = thicknessEdge * (1.0_RKIND - layerCenterSigma(iLevel) )
! Calculate SIA velocity
normalVelocity(iLevel,iEdge) = basalVelocity + &
- 0.5 * ratefactor * (rhoi * gravity)**3 * slopeOnEdge**3 * &
+ 0.5_RKIND * ratefactor * (rhoi * gravity)**3 * slopeOnEdge**3 * &
(thicknessEdge**4 - (thicknessEdge - layerCenterHeightOnEdge)**4)
end do ! Levels
endif
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-04 21:59:39 UTC (rev 2425)
+++ branches/land_ice_projects/implement_core/src/core_land_ice/mpas_land_ice_tendency.F 2013-02-04 22:46:17 UTC (rev 2426)
@@ -119,7 +119,7 @@
nVertLevels = mesh % nVertLevels
! 0 tendency
- layerThickness_tend = 0.0
+ layerThickness_tend = 0.0_RKIND
select case (config_thickness_advection)
case ('FO-Upwind') !===================================================
@@ -176,7 +176,7 @@
! 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
+ mesh % marineBasalMassBal % array = 0.0_RKIND
end where
layerThickness_tend(nVertLevels,:) = layerThickness_tend(nVertLevels,:) &
+ mesh % marineBasalMassBal % array ! (tendency in meters per year)
@@ -257,7 +257,7 @@
integer :: iEdge, k
! 0 tendency
- tracer_tendency = 0.0
+ tracer_tendency = 0.0_RKIND
select case (config_tracer_advection)
case ('FCT') !===================================================
@@ -270,13 +270,13 @@
allocate(uh(mesh % nVertLevels, mesh % nEdges + 1))
! Setup 0 vertical velocity field - vertical velocity not needed with sigma coordinates
- wTop = 0.0
+ wTop = 0.0_RKIND
! prepare for FCT call
! This is a kluge because the FCT code will divide by zero if thickness is zero \todo: deal with more elegantly!
- where (layerThickness == 0.0)
- layerThickness = 1.0e-12
+ where (layerThickness == 0.0_RKIND)
+ layerThickness = 1.0e-12_RKIND
end where
! FCT code needs u*h on edges - this could be calculated elsewhere and saved, potentially (e.g. it's already calculated locally in the FO upwind thickness routine)
@@ -300,8 +300,8 @@
!print *,'tracer_tendency', maxval(tracer_tendency), minval(tracer_tendency)
! Set thickness back to 0 where needed. This is a kluge because the FCT code will divide by zero if thickness is zero \todo: deal with more elegantly!
- where (layerThickness == 1.0e-12)
- layerThickness = 0.0
+ where (layerThickness == 1.0e-12_RKIND)
+ layerThickness = 0.0_RKIND
end where
@@ -370,6 +370,7 @@
lowerSurface, bedTopography, layerThicknessFractions, beta
integer, dimension(:), pointer :: cellMask
real (kind=RKIND), dimension(:,:), pointer :: layerThickness
+ real (kind=RKIND) :: thisThk
nCells = mesh % nCells
nVertLevels = mesh % nVertLevels
@@ -392,7 +393,7 @@
! we get a fresh copy of the array from the annual_forcing subroutine.
! Note: some velocity solvers may do this on their own, but we are doing it here for completeness.
where ( MASK_IS_FLOATING(cellMask) )
- beta = 0.0
+ beta = 0.0_RKIND
end where
! Lower surface is based on floatation for floating ice. For grounded ice it is the bed.
@@ -413,9 +414,12 @@
! set up layerThickness from thickness using the layerThicknessFractions (vertical remapping)
do iCell = 1, nCells
+ thisThk = sum(layerThickness(:,iCell))
do iLevel = 1, nVertLevels
- layerThickness(iLevel,iCell) = layerThicknessFractions(iLevel)*thickness(iCell)
+ layerThickness(iLevel,iCell) = layerThicknessFractions(iLevel) * thisThk
end do
+ ! Check for conservation of mass. Put any residual in the top layer.
+ layerThickness(1,iCell) = layerThickness(1,iCell) + (thisThk - sum(layerThickness(:,iCell)) )
end do
end subroutine land_ice_diagnostic_solve
@@ -494,8 +498,8 @@
cell2 = cellsOnEdge(2,iEdge)
do k=1, nVertLevels
! Calculate h on edges using first order
- VelSign = sign(1.0, normalVelocity(k, iEdge))
- layerThicknessEdge(k,iEdge) = max(VelSign * layerThickness(k, cell1), VelSign * (-1.0) * layerThickness(k, cell2))
+ 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
end do
@@ -708,8 +712,12 @@
real (kind=RKIND) :: flux, VelSign, maxAllowableDt, VelSignSum
real (kind=RKIND), dimension(:), pointer :: dvEdge, areaCell, dcEdge
+ ! Only needed for optional check for mass conservation
+ integer :: iCell
+ real :: tendVolSum
+
err = 0
- MinOfMaxAllowableDt = 1.0e36
+ MinOfMaxAllowableDt = 1.0e36_RKIND
nEdges = grid % nEdges
nVertLevels = grid % nVertLevels
@@ -720,12 +728,8 @@
areaCell => grid % areaCell % array
do iEdge=1,nEdges
- cell1 = cellsOnEdge(1,iEdge)
- cell2 = cellsOnEdge(2,iEdge)
-! if ( (layerThickness(1,cell1) .gt. 0.0) .or. (layerThickness(1,cell2) .gt. 0.0) ) then
- ! Don't calculate for non-ice cells - would result in divide by 0 (just need to check one layer)
- ! \todo: this should use the edgeMask and use the dynamic thickness limit
- if ( MASK_IS_THICK_ICE(edgeMask(iEdge)) ) then
+ cell1 = cellsOnEdge(1,iEdge)
+ cell2 = cellsOnEdge(2,iEdge)
VelSignSum = 0.0
do k=1, nVertLevels
@@ -739,7 +743,7 @@
CellUpwind = cell2
CellDownwind = cell1
endif
- maxAllowableDt = (0.5 * dcEdge(iEdge)) / (abs(normalVelocity(k, iEdge)) + 1.0e-36) ! in years
+ 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
err = err + 1
@@ -747,14 +751,13 @@
MinOfMaxAllowableDt = min(MinOfMaxAllowableDt, maxAllowableDt)
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)
+ 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
endif
- end if
end do
if (err .gt. 0) then
@@ -764,7 +767,14 @@
print *, ' Maximum allowable time step (yrs) on this processor is ~', MinOfMaxAllowableDt
+ ! Optional check for mass conservation
+ !tendVolSum = 0.0_RKIND
+ !do iCell=1, grid % nCells
+ ! tendVolSum = tendVolSum + sum(tend(:,iCell)) * areaCell(iCell)
+ !end do
+ !print *,'SUM OF VOLUME TENDENCY ======', tendVolSum
+
!--------------------------------------------------------------------
end subroutine land_ice_tend_h_layers_fo_upwind!}}}
@@ -861,7 +871,7 @@
! counters, mesh variables, index variables
integer :: iCell, nCells, iEdge, neighborIndex, edgeIndex, k, nVertLevels, jEdge, kEdge
! stuff for making calculations
- real (kind=RKIND) :: numUpwindNeighbors, sumUpwindThickness, numOceanNeighbors, Href, projectedThickness, ExtraVolume, tendencyLayerAdjustment, VelSign, minHref, maxHref
+ real (kind=RKIND) :: numUpwindNeighbors, sumUpwindThickness, numOceanNeighbors, Href, projectedThickness, ExtraVolume, tendencyLayerVolumeRateAdjustment, VelSign, minHref, maxHref
! Only execute the contents of this subroutine if CFBC is turned on.
@@ -947,13 +957,14 @@
! Pass extra mass to neighboring ice-free (ocean only) cells by adjusting tendency in those cells
! I am spreading the extra ice volume evenly amongst all vertical layers. This is a good approximation for SSA flow.
! TODO this could be weighted by the edge length of the various neighbors
- tendencyLayerAdjustment = ExtraVolume / numOceanNeighbors / areaCell(neighborIndex) / nVertLevels / (deltat / SecondsInYear)
+ tendencyLayerVolumeRateAdjustment = ExtraVolume / numOceanNeighbors / nVertLevels / (deltat / SecondsInYear)
! Add the tendency to the neighbor. (Need to add it to the existing value
! in case other overfilled margin cells have already contributed to this neighbor.)
- layerThickness_tend(:,neighborIndex) = layerThickness_tend(:,neighborIndex) + tendencyLayerAdjustment
- ! Lower the tendency in the overfilled cell by the same amount you are passing forward
- layerThickness_tend(:,iCell) = layerThickness_tend(:,iCell) - tendencyLayerAdjustment
+ layerThickness_tend(:,neighborIndex) = layerThickness_tend(:,neighborIndex) + tendencyLayerVolumeRateAdjustment / areaCell(neighborIndex)
+ ! Lower the tendency in the overfilled cell by the same volume you are passing forward
+ layerThickness_tend(:,iCell) = layerThickness_tend(:,iCell) - &
+ tendencyLayerVolumeRateAdjustment / areaCell(iCell)
! === Now calculate a consistent velocity on the edges passing mass to the new cells ===
! Find the edge index (have to loop through edgesOnCell for both cells until we find the matching edge)
@@ -971,9 +982,9 @@
do k = 1, nVertLevels
! These velocities should have been 0 - we are now assigning a velocity
! consistent with the flux we moved (using FO upwinding which is appropriate for advance)
- normalVelocity(k,edgeIndex) = tendencyLayerAdjustment / &
+ normalVelocity(k,edgeIndex) = tendencyLayerVolumeRateAdjustment / &
(dvEdge(edgeIndex) * layerThicknessEdge(k, edgeIndex)) &
- * areaCell(neighborIndex) * VelSign
+ * VelSign
end do
end if
end do
Modified: branches/land_ice_projects/implement_core/src/core_land_ice/mpas_land_ice_time_integration_forwardeuler.F
===================================================================
--- branches/land_ice_projects/implement_core/src/core_land_ice/mpas_land_ice_time_integration_forwardeuler.F 2013-02-04 21:59:39 UTC (rev 2425)
+++ branches/land_ice_projects/implement_core/src/core_land_ice/mpas_land_ice_time_integration_forwardeuler.F 2013-02-04 22:46:17 UTC (rev 2426)
@@ -63,7 +63,6 @@
real (kind=RKIND) :: allowableDt, allowableDtMin
-
! During integration, time level 1 stores the model state at the beginning of the
! time step, and time level 2 stores the state advanced dt in time by timestep(...)
! (time level 1 should not be modified.)
@@ -175,6 +174,7 @@
! Commented out usage for using FCT for thickness
!!!stateNew % sup_thickness % array(1,:,:) = (stateOld % tracers % array(stateOld % index_temperature,:,:) * layerThicknessOld + tracer_tendency(stateOld % index_temperature, :, :) * dt / SecondsInYear) / (layerThicknessNew+1.0e-12)
+
layerThicknessNew = layerThicknessOld + layerThickness_tend * (dt / SecondsInYear)
thicknessNew = sum(layerThicknessNew, 1)
@@ -188,21 +188,21 @@
! if holding advance within initial extent of ice, set thickness to 0 anywhere it has expanded beyond initial extent
if (config_allow_additional_advance .eqv. .false.) then
where ( MASK_WAS_INITIALLY_NOT_ICE(cellMaskOld) )
- thicknessNew = 0.0
+ thicknessNew = 0.0_RKIND
end where
endif
! reset negative thickness to 0. This should not happen unless negative MB is larger than entire ice column.
- where (thicknessNew .lt. 0.0)
+ where (thicknessNew < 0.0_RKIND)
mask = 1
- thicknessNew = 0.0
+ thicknessNew = 0.0_RKIND
end where
- if (sum(mask) .gt. 0) then
+ if (sum(mask) > 0) then
print *,' Cells with negative thickness (set to 0):',sum(mask)
endif
mask = 0
- where (thicknessNew .gt. 0.0)
+ where (thicknessNew > 0.0_RKIND)
mask = 1
end where
print *,' Cells with nonzero thickness:', sum(mask)
@@ -212,12 +212,12 @@
! Calculate new tracer values =================
if (config_tracer_advection .ne. 'None') then
do iTracer = 1, size(tracersNew, 1)
- where (layerThicknessNew .gt. 0.0)
+ where (layerThicknessNew > 0.0_RKIND)
tracersNew(iTracer,:,:) = (tracersOld(iTracer,:,:) * layerThicknessOld &
+ tracer_tendency(iTracer,:,:) * dt / SecondsInYear) / (layerThicknessNew)
elsewhere
! May or may not want to assign tracer values to non-ice cells
- tracersNew(iTracer,:,:) = 0.0
+ tracersNew(iTracer,:,:) = 0.0_RKIND
end where
end do
endif
@@ -227,6 +227,7 @@
! === Calculate diagnostic variables for new state =====================
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
@@ -273,6 +274,7 @@
! Pass in the old velocity to be used as an initial guess for velocity solvers that need one.
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, &
</font>
</pre>