<p><b>mhoffman@lanl.gov</b> 2012-04-18 15:32:11 -0600 (Wed, 18 Apr 2012)</p><p>BRANCH COMMIT<br>
Resolved instability in the thickness advection with the Forward Euler time stepper by (properly) using velocity from previous time step to calculate thickness tendency. Reorganized the Forward Euler time integration subroutine to clearly reflect the proper use of the two time levels. Modified the calls to the velocity solvers to make the velocity variable that is solved to be independent of the state used to calculate velocity.<br>
</p><hr noshade><pre><font color="gray">Modified: branches/land_ice_projects/implement_core/src/core_land_ice/mpas_land_ice_lifev.F
===================================================================
--- branches/land_ice_projects/implement_core/src/core_land_ice/mpas_land_ice_lifev.F 2012-04-18 19:45:23 UTC (rev 1793)
+++ branches/land_ice_projects/implement_core/src/core_land_ice/mpas_land_ice_lifev.F 2012-04-18 21:32:11 UTC (rev 1794)
@@ -289,7 +289,7 @@
!***********************************************************************
- subroutine land_ice_lifev_solve(mesh, state, err)
+ subroutine land_ice_lifev_solve(mesh, state, normalVelocity, err)
!-----------------------------------------------------------------
!
@@ -300,14 +300,19 @@
type (mesh_type), intent(in) :: &
mesh !< Input: mesh information
+ type (state_type), intent(in) :: &
+ state !< Input: state information
+ ! Note: This state is meant to be unmodified (input only).
+
!-----------------------------------------------------------------
!
! input/output variables
!
!-----------------------------------------------------------------
- type (state_type), intent(in) :: &
- state !< Input/Output: state information
+ real (kind=RKIND), intent(inout), dimension(:,:) :: &
+ normalVelocity !< Input/Output: the result of the velocity solve
+ ! Note: this velocity should be used for both input or output - it may be different than the velocity included in state, and the velocity included in state should not be touched. If land_ice_vel_solve() returns anything else besides normalVelocity, those variables should also be treated in this way. This allow the time integration scheme to decide how to treat multiple time levels.
!-----------------------------------------------------------------
!
@@ -324,8 +329,6 @@
!-----------------------------------------------------------------
real (kind=RKIND), dimension(:), pointer :: &
thickness, lowerSurface, beta
- real (kind=RKIND), dimension(:,:), pointer :: &
- normalVelocity
real (kind=RKIND), dimension(:,:,:), pointer :: &
tracers
@@ -337,7 +340,6 @@
lowerSurface => state % lowerSurface % array
beta => state % beta % array
tracers => state % tracers % array
- normalVelocity => state % normalVelocity % array
index_temperature = state % index_temperature
Modified: branches/land_ice_projects/implement_core/src/core_land_ice/mpas_land_ice_sia.F
===================================================================
--- branches/land_ice_projects/implement_core/src/core_land_ice/mpas_land_ice_sia.F 2012-04-18 19:45:23 UTC (rev 1793)
+++ branches/land_ice_projects/implement_core/src/core_land_ice/mpas_land_ice_sia.F 2012-04-18 21:32:11 UTC (rev 1794)
@@ -261,7 +261,7 @@
!***********************************************************************
- subroutine land_ice_sia_solve(mesh, state, err)
+ subroutine land_ice_sia_solve(mesh, state, normalVelocity, err)
use mpas_constants, only: gravity
!-----------------------------------------------------------------
@@ -273,14 +273,19 @@
type (mesh_type), intent(in) :: &
mesh !< Input: mesh information
+ type (state_type), intent(in) :: &
+ state !< Input: state information
+ ! Note: This state is meant to be unmodified (input only).
+
!-----------------------------------------------------------------
!
! input/output variables
!
!-----------------------------------------------------------------
- type (state_type), intent(in) :: &
- state !< Input/Output: state information
+ real (kind=RKIND), intent(inout), dimension(:,:) :: &
+ normalVelocity !< Input/Output: the result of the velocity solve
+ ! Note: this velocity should be used for both input or output - it may be different than the velocity included in state, and the velocity included in state should not be touched. If land_ice_vel_solve() returns anything else besides normalVelocity, those variables should also be treated in this way. This allow the time integration scheme to decide how to treat multiple time levels.
!-----------------------------------------------------------------
!
@@ -300,7 +305,6 @@
real (kind=RKIND), dimension(:), pointer :: dcEdge
real (kind=RKIND), dimension(:,:), pointer :: layerThickness
integer, dimension(:,:), pointer :: cellsOnEdge
- real (kind=RKIND), dimension(:,:), pointer :: normalVelocity
real (kind=RKIND), dimension(mesh%nVertLevels) :: layerFractionAboveBed
integer :: nVertLevels, nCells, nEdges, iCell, iLevel, iEdge
real (kind=RKIND) :: ratefactor, basalVelocity, slopeOnEdge, thicknessOnEdge
@@ -320,7 +324,6 @@
thickness => state % thickness % array
layerThickness => state % layerThickness % array
- normalVelocity => state % normalVelocity % array
rhoi = config_ice_density
@@ -385,6 +388,7 @@
cell1 = cellsOnEdge(1,iEdge)
cell2 = cellsOnEdge(2,iEdge)
! Check if an edge has at least one adjacent cell with ice
+ ! (to turn off advancing margins with FO thickness advection, replace the .or. with .and.
if ( (thickness(cell1).gt.0.0) .or. (thickness(cell2).gt.0.0) ) then
edgeIsIce = .true.
else
Modified: branches/land_ice_projects/implement_core/src/core_land_ice/mpas_land_ice_time_integration.F
===================================================================
--- branches/land_ice_projects/implement_core/src/core_land_ice/mpas_land_ice_time_integration.F 2012-04-18 19:45:23 UTC (rev 1793)
+++ branches/land_ice_projects/implement_core/src/core_land_ice/mpas_land_ice_time_integration.F 2012-04-18 21:32:11 UTC (rev 1794)
@@ -19,7 +19,7 @@
use mpas_configure
use mpas_constants
use mpas_dmpar
- use land_ice_time_integration_ForwardEuler
+ use land_ice_time_integration_forwardeuler
contains
@@ -51,7 +51,7 @@
write(*,*) 'Using ', trim(config_time_integration), ' time integration.'
select case (config_time_integration)
case ('ForwardEuler')
- call land_ice_time_integrator_ForwardEuler(domain, dt)
+ call land_ice_time_integrator_forwardeuler(domain, dt)
case ('RK4')
write(*,*) trim(config_time_integration), ' is not currently supported.'
!call mpas_dmpar_abort(dminfo)
@@ -62,7 +62,7 @@
block => domain % blocklist
do while (associated(block))
- ! Assign the time stamp
+ ! Assign the time stamp for this time step
block % state % time_levs(2) % state % xtime % scalar = timeStamp
! ! Abort the simulation if NaNs occur in the velocity field
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-04-18 19:45:23 UTC (rev 1793)
+++ branches/land_ice_projects/implement_core/src/core_land_ice/mpas_land_ice_time_integration_ForwardEuler.F 2012-04-18 21:32:11 UTC (rev 1794)
@@ -1,6 +1,6 @@
!|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
!
-! land_ice_time_integration_ForwardEuler
+! land_ice_time_integration_forwardeuler
!
!> \brief MPAS land ice Forward Euler time integration schem
!> \author Matt Hoffman
@@ -12,7 +12,7 @@
!-----------------------------------------------------------------------
-module land_ice_time_integration_ForwardEuler
+module land_ice_time_integration_forwardeuler
use mpas_vector_reconstruction
use mpas_grid_types
@@ -24,7 +24,7 @@
contains
- subroutine land_ice_time_integrator_ForwardEuler(domain, dt)
+ subroutine land_ice_time_integrator_forwardeuler(domain, dt)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Advance model state forward in time by the specified time step using
! Forward Euler method
@@ -41,11 +41,11 @@
real (kind=RKIND), intent(in) :: dt
type (block_type), pointer :: block
- type (state_type), pointer :: state
+ type (state_type), pointer :: stateOld, stateNew
type (mesh_type), pointer :: mesh
- real (kind=RKIND), dimension(:), pointer :: thickness, thickness_tend, layerThicknessFractions
- real (kind=RKIND), dimension(:,:), pointer :: normalVelocity, layerThickness !, layerThickness_tend
+ real (kind=RKIND), dimension(:), pointer :: thicknessOld, thicknessNew, thickness_tend, layerThicknessFractions
+ real (kind=RKIND), dimension(:,:), pointer :: normalVelocityOld, normalVelocityNew, layerThicknessOld, layerThicknessNew !, layerThickness_tend
real (kind=RKIND), dimension(:,:), pointer :: uReconstructX, uReconstructY, uReconstructZ, &
uReconstructZonal, uReconstructMeridional
@@ -55,72 +55,105 @@
integer :: err
+ ! 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.)
block => domain % blocklist
do while (associated(block))
+! === Initialize pointers ========================
+ ! Mesh information
mesh => block % mesh
nVertLevels = mesh % nVertLevels
layerThicknessFractions => mesh % layerThicknessFractions % array
- state => block % state % time_levs(2) % state
- thickness => state % thickness % array
- layerThickness => state % layerThickness % array
- normalVelocity => state % normalVelocity % array
+ ! State at time n
+ stateOld => block % state % time_levs(1) % state
+ thicknessOld => stateOld % thickness % array
+ layerThicknessOld => stateOld % layerThickness % array
+ normalVelocityOld => stateOld % normalVelocity % array
- uReconstructX => state % uReconstructX % array
- uReconstructY => state % uReconstructY % array
- uReconstructZ => state % uReconstructZ % array
- uReconstructZonal => state % uReconstructZonal % array
- uReconstructMeridional => state % uReconstructMeridional % 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
+ uReconstructX => stateNew % uReconstructX % array
+ uReconstructY => stateNew % uReconstructY % array
+ uReconstructZ => stateNew % uReconstructZ % array
+ uReconstructZonal => stateNew % uReconstructZonal % array
+ uReconstructMeridional => stateNew % uReconstructMeridional % array
+ ! Tendencies
!layerThickness_tend => block % tend % layerThickness % array
thickness_tend => block % tend % thickness % array
- k = size(thickness)
- allocate(mask(k))
+ allocate(mask(mesh%nCells + 1))
- ! Solve the velocity
- call land_ice_vel_solve(mesh, block % state % time_levs(1) % state, err)
- block % state % time_levs(2) % state % normalVelocity % array = block % state % time_levs(1) % state % normalVelocity % array
- !state % xtime % scalar = timeStamp
-
- ! Update the thickness and tracers (move to another module eventually)
- call lice_tend_h(mesh, normalVelocity, layerThickness, thickness, thickness_tend, dt, err)
- thickness = thickness + thickness_tend * dt / SecondsInYear
- ! Commented lines below calculate the thickness tendency per layer instead of for the whole column
+! === Calculate Tendencies ========================
+ ! Calculate thickness tendency using state at time n
+ call lice_tend_h(mesh, normalVelocityOld, layerThicknessOld, thicknessOld, thickness_tend, dt, err)
+
+ ! update halos on thickness tend
+ !call ()
+
+ ! (Calculate tracer tendencies)
+ !call ()
+
+
+ ! Commented lines below calculate the thickness tendency per layer instead of for the whole column (last 2 lines should be moved to following 'section' of code)
!call lice_tend_h_layers(mesh, normalVelocity, layerThickness, layerThickness_tend, dt, err)
!layerThickness = layerThickness + layerThickness_tend * dt / SecondsInYear
!thickness = sum(layerThickness, 1)
- ! Remap layer thicknesses
- do k = 1, nVertLevels
- layerThickness(k, :) = thickness * layerThicknessFractions(k)
- end do
+! === Compute new state for prognostic variables ==================================
+
+ ! Update thickness (and tracers eventually)
+ thicknessNew = thicknessOld + thickness_tend * dt / SecondsInYear
+
!Optionally print some information about the new state
- print *, 'thickness maxval:', maxval(thickness)
+ print *, 'thickness maxval:', maxval(thicknessNew)
mask = 0
- where (thickness .lt. 0.0)
+ where (thicknessNew .lt. 0.0)
mask = 1
- thickness = 0.0
+ thicknessNew = 0.0
end where
- print *,'Cells with negative thickness (set to 0):',sum(mask)
+ if (sum(mask) .gt. 0) then
+ print *,'Cells with negative thickness (set to 0):',sum(mask)
+ endif
mask = 0
- where (thickness .gt. 0.0)
+ where (thicknessNew .gt. 0.0)
mask = 1
end where
print *,'Cells with nonzero thickness:', sum(mask)
+ ! Solve the velocity (using old state)
+ call land_ice_vel_solve(mesh, stateOld, normalVelocityNew, err)
+ ! update halos on velocity
+ ! call ()
+
+! === Calculate diagnostic variables for new state =====================
+ ! Remap layer thicknesses (could be moved to diagnostic solve)
+ do k = 1, nVertLevels
+ layerThicknessNew(k, :) = thicknessNew * layerThicknessFractions(k)
+ end do
+
! call land_ice_diagnostic_solve()
- ! Reconstruct velocities for this time step
- call mpas_reconstruct(mesh, normalVelocity, &
+ ! Reconstruct velocities for this time step (could be moved to diagnostic solve)
+ call mpas_reconstruct(mesh, normalVelocityNew, &
uReconstructX, uReconstructY, uReconstructZ, &
uReconstructZonal, uReconstructMeridional)
+! === Cleanup & Misc. =============================
+
+
+! ================================
+
block => block % next
end do
@@ -130,6 +163,6 @@
-end module land_ice_time_integration_ForwardEuler
+end module land_ice_time_integration_forwardeuler
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 2012-04-18 19:45:23 UTC (rev 1793)
+++ branches/land_ice_projects/implement_core/src/core_land_ice/mpas_land_ice_vel.F 2012-04-18 21:32:11 UTC (rev 1794)
@@ -238,7 +238,7 @@
!***********************************************************************
- subroutine land_ice_vel_solve(mesh, state, err)
+ subroutine land_ice_vel_solve(mesh, state, normalVelocity, err)
use land_ice_SIA
@@ -251,14 +251,19 @@
type (mesh_type), intent(in) :: &
mesh !< Input: mesh information
+ type (state_type), intent(in) :: &
+ state !< Input: state information
+ ! Note: This state is meant to be unmodified (input only).
+
!-----------------------------------------------------------------
!
! input/output variables
!
!-----------------------------------------------------------------
- type (state_type), intent(in) :: &
- state !< Input/Output: state information
+ real (kind=RKIND), intent(inout), dimension(:,:) :: &
+ normalVelocity !< Input/Output: the result of the velocity solve
+ ! Note: this velocity should be used for both input or output - it may be different than the velocity included in state, and the velocity included in state should not be touched. If land_ice_vel_solve() returns anything else besides normalVelocity, those variables should also be treated in this way. This allow the time integration scheme to decide how to treat multiple time levels.
!-----------------------------------------------------------------
!
@@ -278,9 +283,9 @@
select case (config_dycore)
case ('L1L2')
- call land_ice_lifev_solve(mesh, state, err)
+ call land_ice_lifev_solve(mesh, state, normalVelocity, err)
case ('SIA')
- call land_ice_sia_solve(mesh, state, err)
+ call land_ice_sia_solve(mesh, state, normalVelocity, err)
case default
write(*,*) trim(config_dycore), ' is not a valid land ice dycore option.'
err = 1
</font>
</pre>