<p><b>mhoffman@lanl.gov</b> 2012-06-01 14:17:58 -0600 (Fri, 01 Jun 2012)</p><p> BRANCH COMMIT -- land ice<br>
Adding SMB back into the thickness tendency calculation. This formulation may break tracer advection (not tested).<br>
Also made the Forward Euler timestepper abort on a CFL violation and suggest an appropriate time step.<br>
</p><hr noshade><pre><font color="gray">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-06-01 19:50:54 UTC (rev 1954)
+++ branches/land_ice_projects/implement_core/src/core_land_ice/mpas_land_ice_tendency.F 2012-06-01 20:17:58 UTC (rev 1955)
@@ -147,10 +147,20 @@
end select !===================================================
- ! Add surface mass balance to tendency
- ! NOTE: 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?)
- !thickness_tend = thickness_tend + mesh % sfcMassBal % array * dt / SecondsInYear
+ ! Add the MB to the tendencies
+ select case (config_thickness_advection)
+ case ('None') !===================================================
+ ! Do nothing - don't add the MB
+ case default
+ ! Add surface mass balance to tendency
+ ! NOTE: 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)
+ ! THIS MIGHT RESULT IN NEGATIVE LAYER THICKNESS!
+ ! Do basal mass balance too?
+
+ end select
+
end subroutine land_ice_tendency_h
@@ -473,10 +483,11 @@
integer, dimension(:,:), pointer :: cellsOnEdge
- real (kind=RKIND) :: flux, VelSign, h_edge, ubar
+ real (kind=RKIND) :: flux, VelSign, h_edge, ubar, maxAllowableDt
real (kind=RKIND), dimension(:), pointer :: dvEdge, areaCell, dcEdge
err = 0
+ maxAllowableDt = 1.0e36
nEdges = grid % nEdges
nVertLevels = grid % nVertLevels
@@ -511,13 +522,19 @@
enddo
ubar = flux / thickness(cellUpwind)
if ( (abs(ubar) * dt/SecondsInYear) .gt. (0.5 * dcEdge(iEdge))) then
+ !maxAllowableDt = min(maxAllowableDt, (0.5 * dcEdge)/ubar )
write(0,*) 'CFL violation at edge', iEdge
- !err = 1
+ err = err + 1
endif
thickness_tend(cellUpwind) = thickness_tend(cellUpwind) - flux * dvEdge(iEdge) / areaCell(cellUpwind)
thickness_tend(cellDownwind) = thickness_tend(cellDownwind) + flux * dvEdge(iEdge) / areaCell(cellDownwind)
endif
end do
+
+ if (err .gt. 0) then
+ write(6,*) 'CFL violation at ', err, ' edges! Maximum time step should be ', maxAllowableDt
+ err = 1
+ endif
!--------------------------------------------------------------------
end subroutine land_ice_tend_h_fo_upwind !}}}
@@ -587,10 +604,11 @@
integer, dimension(:,:), pointer :: cellsOnEdge
- real (kind=RKIND) :: flux, VelSign
+ real (kind=RKIND) :: flux, VelSign, maxAllowableDt, MinOfMaxAllowableDt
real (kind=RKIND), dimension(:), pointer :: dvEdge, areaCell, dcEdge
err = 0
+ MinOfMaxAllowableDt = 1.0e36
nEdges = grid % nEdges
nVertLevels = grid % nVertLevels
@@ -613,9 +631,13 @@
endif
if (layerThickness(1,cellUpwind) .gt. 0.0) then ! Don't calculate for non-ice cells - would result in divide by 0 (just need to check one layer)
do k=1, nVertLevels
- if ( (abs(normalVelocity(k, iEdge)) * dt/SecondsInYear) .gt. (0.5 * dcEdge(iEdge))) then
- print *, 'CFL violation at level, edge', k, iEdge
+ maxAllowableDt = (0.5 * dcEdge(iEdge)) / (abs(normalVelocity(k, iEdge)) + 1.0e-36) ! in years
+ if ( (dt/SecondsInYear) .gt. maxAllowableDt ) then
+ !write(0,*) 'CFL violation at level, edge', k, iEdge
+ err = err + 1
endif
+ 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)
@@ -623,6 +645,13 @@
end if
end do
+ if (err .gt. 0) then
+ write(0,*) 'CFL violation on ', err, ' level-edges! Maximum time step (yrs) should be ', MinOfMaxAllowableDt
+ err = 1
+ endif
+
+ print *, ' Maximum allowable time step (yrs) is ~', MinOfMaxAllowableDt
+
!--------------------------------------------------------------------
end subroutine land_ice_tend_h_layers_fo_upwind!}}}
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-06-01 19:50:54 UTC (rev 1954)
+++ branches/land_ice_projects/implement_core/src/core_land_ice/mpas_land_ice_time_integration_forwardeuler.F 2012-06-01 20:17:58 UTC (rev 1955)
@@ -69,6 +69,9 @@
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
@@ -114,12 +117,23 @@
! 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
+
! update halos on tracer tend (perhaps move to land_ice_tendency_tracers after we've updated the halo update calls)
- call mpas_dmpar_exch_halo_field3d_real(dminfo, tracer_tendency, &
+ 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
+
call mpas_timer_stop("calculate tendencies")
@@ -134,7 +148,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
+ layerThicknessNew = layerThicknessOld + layerThickness_tend * (dt / SecondsInYear)
thicknessNew = sum(layerThicknessNew, 1)
!Optionally print some information about the new thickness
@@ -142,6 +156,8 @@
!print *, 'thicknessOld maxval:', maxval(thicknessOld(1:mesh % nCellsSolve))
print *, ' thicknessNew maxval:', maxval(thicknessNew(1:mesh % nCellsSolve))
mask = 0
+
+ ! reset negative thickness to 0. This should not happen unless negative MB is larger than entire ice column.
where (thicknessNew .lt. 0.0)
mask = 1
thicknessNew = 0.0
</font>
</pre>