<p><b>mhoffman@lanl.gov</b> 2012-04-17 14:15:32 -0600 (Tue, 17 Apr 2012)</p><p>BRANCH COMMIT<br>
Added basic time-stepping (Forward Euler) and thickness advection (FO Upwinding). Works for dome test case, but thickness instabilities appear after about 10 years.<br>
Also cleaned up Registry and default namelist files. Changed the way dt is specified in the namelist file to allow it to be specified in years (rather than seconds).<br>
</p><hr noshade><pre><font color="gray">Modified: branches/land_ice_projects/implement_core/namelist.input.land_ice
===================================================================
--- branches/land_ice_projects/implement_core/namelist.input.land_ice 2012-04-17 18:59:11 UTC (rev 1789)
+++ branches/land_ice_projects/implement_core/namelist.input.land_ice 2012-04-17 20:15:32 UTC (rev 1790)
@@ -1,16 +1,9 @@
&land_ice_model
- config_test_case = 1
- config_time_integration = 'RK4'
- config_dt = 60.0
+ config_time_integration = 'ForwardEuler'
+ config_dycore = 'SIA'
+ config_dt_years = 1.0
config_start_time = '0000-01-01_00:00:00'
- config_stop_time = '0000-01-01_00:01:00'
- config_run_duration = '00:01:00'
- config_stats_interval = 0
- config_h_ScaleWithMesh = .false.
- config_thickness_adv_order = 2
- config_tracer_adv_order = 2
- config_positive_definite = .false.
- config_monotonic = .false.
+ config_stop_time = '0001-01-01_00:00:00'
/
&io
Modified: branches/land_ice_projects/implement_core/src/core_land_ice/Makefile
===================================================================
--- branches/land_ice_projects/implement_core/src/core_land_ice/Makefile 2012-04-17 18:59:11 UTC (rev 1789)
+++ branches/land_ice_projects/implement_core/src/core_land_ice/Makefile 2012-04-17 20:15:32 UTC (rev 1790)
@@ -1,7 +1,9 @@
.SUFFIXES: .F .o
OBJS = mpas_land_ice_mpas_core.o \
+ mpas_land_ice_tendency.o \
mpas_land_ice_time_integration.o \
+ mpas_land_ice_time_integration_ForwardEuler.o \
mpas_land_ice_vel.o \
mpas_land_ice_lifev.o \
mpas_land_ice_sia.o \
@@ -12,17 +14,21 @@
core_land_ice: $(OBJS)
ar -ru libdycore.a $(OBJS)
-mpas_land_ice_time_integration.o: mpas_land_ice_vel.o
+mpas_land_ice_time_integration.o: mpas_land_ice_time_integration_ForwardEuler.o
+mpas_land_ice_time_integration_ForwardEuler.o: mpas_land_ice_vel.o mpas_land_ice_tendency.o
+
mpas_land_ice_vel.o: mpas_land_ice_lifev.o mpas_land_ice_sia.o
+mpas_land_ice_tendency.o:
+
mpas_land_ice_lifev.o:
mpas_land_ice_sia.o:
mpas_land_ice_global_diagnostics.o:
-mpas_land_ice_mpas_core.o: mpas_land_ice_global_diagnostics.o mpas_land_ice_time_integration.o mpas_land_ice_vel.o
+mpas_land_ice_mpas_core.o: mpas_land_ice_global_diagnostics.o mpas_land_ice_time_integration.o
clean:
$(RM) *.o *.mod *.f90 libdycore.a
Modified: branches/land_ice_projects/implement_core/src/core_land_ice/Registry
===================================================================
--- branches/land_ice_projects/implement_core/src/core_land_ice/Registry 2012-04-17 18:59:11 UTC (rev 1789)
+++ branches/land_ice_projects/implement_core/src/core_land_ice/Registry 2012-04-17 20:15:32 UTC (rev 1790)
@@ -3,21 +3,28 @@
%
% namelist type namelist_record name default_value
%
-namelist integer land_ice_model config_test_case 5
-namelist character land_ice_model config_time_integration RK4
-namelist real land_ice_model config_dt 172.8
-namelist character land_ice_model config_calendar_type 360day
+% Land ice core will use dt in years if it is supplied. Otherwise it will use the dt in seconds. Default config_dt_seconds = 1 yr.
+namelist real land_ice_model config_dt_years 0.0
+namelist real land_ice_model config_dt_seconds 31536000.0
+namelist character land_ice_model config_calendar_type gregorian_noleap
+% Valid calendar types are: gregorian_noleap, gregorian, and 360day
namelist character land_ice_model config_start_time 0000-01-01_00:00:00
namelist character land_ice_model config_stop_time none
namelist character land_ice_model config_run_duration none
+namelist character land_ice_model config_dycore SIA
+% valid dycores are SIA, L1L2, HO, Stokes. All but SIA require compiling with LifeV libraries.
+namelist character land_ice_model config_time_integration ForwardEuler
+% valid time integration is ForwardEuler
+% %namelist character land_ice_model config_advection FO-Upwind
+%% valid time integration is FO-Upwind - this could be added once other advection schemes are added.
+namelist real land_ice_model config_ice_density 900.0
+% The following options are used by the ocean core and may be useful in the future
namelist integer land_ice_model config_stats_interval 100
namelist logical land_ice_model config_h_ScaleWithMesh false
-namelist integer land_ice_model config_thickness_adv_order 2
-namelist integer land_ice_model config_tracer_adv_order 2
-namelist logical land_ice_model config_positive_definite false
-namelist logical land_ice_model config_monotonic false
-namelist character land_ice_model config_dycore SIA
-namelist real land_ice_model config_ice_density 900.0
+%namelist integer land_ice_model config_thickness_adv_order 1
+%namelist integer land_ice_model config_tracer_adv_order 2
+%namelist logical land_ice_model config_positive_definite false
+%namelist logical land_ice_model config_monotonic false
namelist character io config_input_name grid.nc
namelist character io config_output_name output.nc
namelist character io config_restart_name restart.nc
@@ -28,7 +35,7 @@
namelist character restart config_restart_interval none
%
-% dim type name_in_file name_in_code
+% dim name_in_file name_in_code
%
dim nCells nCells
dim nEdges nEdges
@@ -122,6 +129,7 @@
% Tendency variables
var persistent real tend_layerThickness ( nVertLevels nCells Time ) 1 - layerThickness tend - -
+var persistent real tend_thickness ( nCells Time ) 1 - thickness tend - -
var persistent real tend_temperature ( nVertLevels nCells Time ) 1 - temperature tend tracers dynamics
% Diagnostic fields: only written to output
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-17 18:59:11 UTC (rev 1789)
+++ branches/land_ice_projects/implement_core/src/core_land_ice/mpas_land_ice_lifev.F 2012-04-17 20:15:32 UTC (rev 1790)
@@ -9,6 +9,7 @@
!> \details
!> This module contains the routines for interfacing with LifeV
!>
+!> MJH: added interfaces following http://people.sc.fsu.edu/~jburkardt/f_src/f90_calls_c++/f90_calls_c++.html
!
!-----------------------------------------------------------------------
@@ -17,6 +18,7 @@
use mpas_grid_types
use mpas_configure
use mpas_dmpar
+ use, intrinsic :: iso_c_binding
implicit none
private
@@ -110,6 +112,21 @@
logical :: keepVertex
+ !-----------------------------------------------------------------
+ !
+ ! LifeV C++ interfaces
+ !
+ !-----------------------------------------------------------------
+
+ interface
+ subroutine velocity_solver_init_l1l2(nVertLevels) bind ( c )
+ use iso_c_binding
+ integer ( c_int ), VALUE :: nVertLevels
+ end subroutine velocity_solver_init_l1l2
+ end interface
+
+
+
err = 0
! Note: the following code assumes that blocklist contains a single block
@@ -158,7 +175,7 @@
!deallocate(verticesMask)
- !call lifev_initialize_L1L2_solver(nVertLevels)
+ !call velocity_solver_init_l1l2(nVertLevels) ! LifeV C++ call
!--------------------------------------------------------------------
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 2012-04-17 18:59:11 UTC (rev 1789)
+++ branches/land_ice_projects/implement_core/src/core_land_ice/mpas_land_ice_mpas_core.F 2012-04-17 20:15:32 UTC (rev 1790)
@@ -12,12 +12,15 @@
integer, parameter :: restartAlarmID = 2
!integer, parameter :: statsAlarmID = 3
+ real(kind=RKIND) :: dt
+
contains
subroutine mpas_core_init(domain, startTimeStamp)
use mpas_configure
use mpas_grid_types
+ use mpas_constants
use land_ice_vel, only: land_ice_vel_init
implicit none
@@ -25,7 +28,6 @@
type (domain_type), intent(inout) :: domain
character(len=*), intent(out) :: startTimeStamp
- real (kind=RKIND) :: dt
type (block_type), pointer :: block
integer :: i, err
@@ -36,7 +38,12 @@
!
! Initialize core
!
- dt = config_dt
+ ! If dt in years is supplied, use it. Otherwise use dt in seconds
+ if (config_dt_years .eq. 0.0) then
+ dt = config_dt_seconds
+ else
+ dt = config_dt_years * SecondsInYear
+ endif
call simulation_clock_init(domain, dt, startTimeStamp)
@@ -212,15 +219,12 @@
integer, intent(inout) :: output_frame
integer :: itimestep
- real (kind=RKIND) :: dt
type (block_type), pointer :: block_ptr
type (MPAS_Time_Type) :: currTime
character(len=32) :: timeStamp
integer :: ierr
- ! Eventually, dt should be domain specific
- dt = config_dt
currTime = mpas_get_clock_time(clock, MPAS_NOW, ierr)
call mpas_get_time(curr_time=currTime, dateTimeString=timeStamp, ierr=ierr)
Added: 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 (rev 0)
+++ branches/land_ice_projects/implement_core/src/core_land_ice/mpas_land_ice_tendency.F 2012-04-17 20:15:32 UTC (rev 1790)
@@ -0,0 +1,269 @@
+!|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
+!
+! land_ice_tendency
+!
+!> \brief MPAS land ice tendency driver
+!> \author Matt Hoffman
+!> \date 17 April 2011
+!> \version SVN:$Id:$
+!> \details
+!> This module contains the routines for computing
+!> various tendencies for land ice. As well as routines
+!> for computing diagnostic variables.
+!
+!-----------------------------------------------------------------------
+
+module land_ice_tendency
+
+ use mpas_grid_types
+ use mpas_configure
+ use mpas_constants
+ use mpas_dmpar
+
+ implicit none
+ private
+ save
+
+ !--------------------------------------------------------------------
+ !
+ ! Public parameters
+ !
+ !--------------------------------------------------------------------
+
+ !--------------------------------------------------------------------
+ !
+ ! Public member functions
+ !
+ !--------------------------------------------------------------------
+ public :: lice_tend_h, &
+ lice_tend_h_layers
+
+ !--------------------------------------------------------------------
+ !
+ ! Private module variables
+ !
+ !--------------------------------------------------------------------
+
+
+
+ contains
+
+
+
+!***********************************************************************
+!
+! subroutine lice_tend_h
+!
+!> \brief Computes tendency term from horizontal advection of thickness
+!> \author Matt Hoffman
+!> \date 16 April 2012
+!> \version SVN:$Id$
+!> \details
+!> This routine computes the horizontal advection tendency for
+!> thickness based on current state and user choices of forcings. Based on
+!> ocn_thick_hadv_tend in the ocean core.
+!
+!-----------------------------------------------------------------------
+
+ subroutine lice_tend_h(grid, normalVelocity, layerThickness, thickness, thickness_tend, dt, err)!{{{
+
+ !-----------------------------------------------------------------
+ !
+ ! input variables
+ !
+ !-----------------------------------------------------------------
+
+ real (kind=RKIND), dimension(:,:), intent(in) :: &
+ normalVelocity !< Input: velocity
+
+ real (kind=RKIND), dimension(:,:), intent(in) :: &
+ layerThickness !< Input: thickness of each layer
+
+ real (kind=RKIND), dimension(:), intent(in) :: &
+ thickness !< Input: thickness
+
+ type (mesh_type), intent(in) :: &
+ grid !< Input: grid information
+
+ real (kind=RKIND) :: dt
+
+ !-----------------------------------------------------------------
+ !
+ ! input/output variables
+ !
+ !-----------------------------------------------------------------
+
+ real (kind=RKIND), dimension(:), intent(inout) :: &
+ thickness_tend !< Input/Output: thickness tendency
+
+ !-----------------------------------------------------------------
+ !
+ ! output variables
+ !
+ !-----------------------------------------------------------------
+
+ integer, intent(out) :: err !< Output: error flag
+
+ !-----------------------------------------------------------------
+ !
+ ! local variables
+ !
+ !-----------------------------------------------------------------
+
+ integer :: iEdge, nEdges, cell1, cell2, nVertLevels, k, CellUpwind, CellDownwind
+
+ integer, dimension(:,:), pointer :: cellsOnEdge
+
+ real (kind=RKIND) :: flux, VelSign, h_edge, ubar
+ real (kind=RKIND), dimension(:), pointer :: dvEdge, areaCell, dcEdge
+
+ err = 0
+
+ nEdges = grid % nEdges
+ nVertLevels = grid % nVertLevels
+ cellsOnEdge => grid % cellsOnEdge % array
+ dvEdge => grid % dvEdge % array
+ dcEdge => grid % dcEdge % array
+ areaCell => grid % areaCell % array
+
+ ! Zero the tendency before accumulating
+ thickness_tend = 0.0
+
+ do iEdge=1,nEdges
+ cell1 = cellsOnEdge(1,iEdge)
+ cell2 = cellsOnEdge(2,iEdge)
+
+ VelSign = sign(1.0, normalVelocity(1, iEdge))
+ if (VelSign .gt. 0.0) then
+ CellUpwind = cell1
+ CellDownwind = cell2
+ else
+ CellUpwind = cell2
+ CellDownwind = cell1
+ endif
+ if (thickness(cellUpwind) .gt. 0.0) then ! Don't calculate for non-ice cells - would result in divide by 0
+ h_edge = thickness(CellUpwind)
+ flux = 0.0
+ do k=1, nVertLevels
+ ! Calculate thickness averaged velocity - this make be calculated externally and passed in
+ flux = flux + layerThickness(k, cellUpwind) * abs(normalVelocity(k, iEdge))
+ enddo
+ ubar = flux / thickness(cellUpwind)
+ if ( (abs(ubar) * dt/SecondsInYear) .gt. (0.5 * dcEdge(iEdge))) then
+ print *, 'CFL violation at edge', iEdge
+ !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
+ !--------------------------------------------------------------------
+
+ end subroutine lice_tend_h!}}}
+
+
+
+!***********************************************************************
+!
+! subroutine lice_tend_h_layers
+!
+!> \brief Computes tendency term from horizontal advection of thickness layers
+!> \author Matt Hoffman
+!> \date 16 April 2012
+!> \version SVN:$Id$
+!> \details
+!> This routine computes the horizontal advection tendency for each
+!> thickness layer based on current state and user choices of forcings. Based on
+!> ocn_thick_hadv_tend in the ocean core. This is an alternative to lice_tend_h
+!> that caclculates the tendency for each layer, which would then need to be
+!> added up to calculate the change in thickness. The two methods yield identical
+!> results. Currently, this is not used but I'm leaving it here in case it is useful later.
+!
+!-----------------------------------------------------------------------
+ subroutine lice_tend_h_layers(grid, normalVelocity, layerThickness, tend, dt, err)!{{{
+
+ !-----------------------------------------------------------------
+ !
+ ! input variables
+ !
+ !-----------------------------------------------------------------
+
+ real (kind=RKIND), dimension(:,:), intent(in) :: &
+ normalVelocity !< Input: velocity
+
+ real (kind=RKIND), dimension(:,:), intent(in) :: &
+ layerThickness !< Input: thickness of each layer at edge
+
+ type (mesh_type), intent(in) :: &
+ grid !< Input: grid information
+
+ real (kind=RKIND) :: dt
+
+ !-----------------------------------------------------------------
+ !
+ ! input/output variables
+ !
+ !-----------------------------------------------------------------
+
+ real (kind=RKIND), dimension(:,:), intent(inout) :: &
+ tend !< Input/Output: velocity tendency
+
+ !-----------------------------------------------------------------
+ !
+ ! output variables
+ !
+ !-----------------------------------------------------------------
+
+ integer, intent(out) :: err !< Output: error flag
+
+ !-----------------------------------------------------------------
+ !
+ ! local variables
+ !
+ !-----------------------------------------------------------------
+
+ integer :: iEdge, nEdges, cell1, cell2, nVertLevels, k
+
+ integer, dimension(:,:), pointer :: cellsOnEdge
+
+ real (kind=RKIND) :: flux, VelSign, h_edge
+ real (kind=RKIND), dimension(:), pointer :: dvEdge, areaCell, dcEdge
+
+ err = 0
+
+ nEdges = grid % nEdges
+ nVertLevels = grid % nVertLevels
+
+ cellsOnEdge => grid % cellsOnEdge % array
+ dvEdge => grid % dvEdge % array
+ dcEdge => grid % dcEdge % array
+ areaCell => grid % areaCell % array
+
+ ! Zero the tendency before accumulating
+ tend = 0.0
+
+ do iEdge=1,nEdges
+ cell1 = cellsOnEdge(1,iEdge)
+ cell2 = cellsOnEdge(2,iEdge)
+ 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
+ endif
+ ! Calculate h on edges using first order
+ VelSign = sign(1.0, normalVelocity(k, iEdge))
+ h_edge = max(0.0, VelSign) * layerThickness(k, cell1) &
+ - min(0.0, VelSign) * layerThickness(k, cell2)
+ !h_edge = (layerThickness(k, cell1) +layerThickness(k, cell2) ) / 2.0
+ flux = normalVelocity(k,iEdge) * dvEdge(iEdge) * h_edge
+ tend(k,cell1) = tend(k,cell1) - flux / areaCell(cell1)
+ tend(k,cell2) = tend(k,cell2) + flux / areaCell(cell2)
+ end do
+ end do
+
+ !--------------------------------------------------------------------
+
+ end subroutine lice_tend_h_layers!}}}
+
+
+
+end module land_ice_tendency
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-17 18:59:11 UTC (rev 1789)
+++ branches/land_ice_projects/implement_core/src/core_land_ice/mpas_land_ice_time_integration.F 2012-04-17 20:15:32 UTC (rev 1790)
@@ -1,3 +1,17 @@
+!|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
+!
+! land_ice_time_integration
+!
+!> \brief MPAS land ice time integration driver
+!> \author Matt Hoffman
+!> \date 17 April 2011
+!> \version SVN:$Id:$
+!> \details
+!> This module contains the main driver routine for calling
+!> the time integration scheme
+!
+!-----------------------------------------------------------------------
+
module land_ice_time_integration
use mpas_vector_reconstruction
@@ -5,7 +19,7 @@
use mpas_configure
use mpas_constants
use mpas_dmpar
- use land_ice_vel, only: land_ice_vel_solve
+ use land_ice_time_integration_ForwardEuler
contains
@@ -34,34 +48,29 @@
integer :: err
-! if (trim(config_time_integration) == 'RK4') then
-! call land_ice_rk4(domain, dt)
-! else
-! write(0,*) 'Unknown time integration option '//trim(config_time_integration)
-! write(0,*) 'Currently, only ''RK4'' is supported.'
-! stop
-! end if
+ write(*,*) 'Using ', trim(config_time_integration), ' time integration.'
+ select case (config_time_integration)
+ case ('ForwardEuler')
+ call land_ice_time_integrator_ForwardEuler(domain, dt)
+ case ('RK4')
+ write(*,*) trim(config_time_integration), ' is not currently supported.'
+ !call mpas_dmpar_abort(dminfo)
+ case default
+ write(*,*) trim(config_time_integration), ' is not a valid land ice time integration option.'
+ !call mpas_dmpar_abort(dminfo)
+ end select
- ! a temporary test that calls the velocity solver on the state at the current time
-
block => domain % blocklist
do while (associated(block))
- state => block % state % time_levs(2) % state
- normalVelocity => block % state % time_levs(2) % state % normalVelocity % array
- uReconstructX => block % state % time_levs(2) % state % uReconstructX % array
- uReconstructY => block % state % time_levs(2) % state % uReconstructY % array
- uReconstructZ => block % state % time_levs(2) % state % uReconstructZ % array
- uReconstructZonal => block % state % time_levs(2) % state % uReconstructZonal % array
- uReconstructMeridional => block % state % time_levs(2) % state % uReconstructMeridional % array
+ ! Assign the time stamp
+ block % state % time_levs(2) % state % xtime % scalar = timeStamp
- call land_ice_vel_solve(block % mesh, state, err)
- state % xtime % scalar = timeStamp
+! ! Abort the simulation if NaNs occur in the velocity field
+! if (isNaN(sum(block % state % time_levs(2) % state % u % array))) then
+! write(0,*) 'Abort: NaN detected'
+! call mpas_dmpar_abort(dminfo)
+! endif
- ! This call might make more sense in mpas_core_run but there already is a block loop here.
- call mpas_reconstruct(block % mesh, normalVelocity, &
- uReconstructX, uReconstructY, uReconstructZ, &
- uReconstructZonal, uReconstructMeridional)
-
block => block % next
end do
Added: 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 (rev 0)
+++ branches/land_ice_projects/implement_core/src/core_land_ice/mpas_land_ice_time_integration_ForwardEuler.F 2012-04-17 20:15:32 UTC (rev 1790)
@@ -0,0 +1,135 @@
+!|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
+!
+! land_ice_time_integration_ForwardEuler
+!
+!> \brief MPAS land ice Forward Euler time integration schem
+!> \author Matt Hoffman
+!> \date 17 April 2011
+!> \version SVN:$Id:$
+!> \details
+!> This module contains the Forware Euler time integration scheme
+!
+!-----------------------------------------------------------------------
+
+
+module land_ice_time_integration_ForwardEuler
+
+ use mpas_vector_reconstruction
+ use mpas_grid_types
+ use mpas_configure
+ use mpas_constants
+ use mpas_dmpar
+ use land_ice_vel, only: land_ice_vel_solve
+ use land_ice_tendency
+
+ contains
+
+ subroutine land_ice_time_integrator_ForwardEuler(domain, dt)
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ ! Advance model state forward in time by the specified time step using
+ ! Forward Euler method
+ !
+ ! Input: domain - current model state in time level 1 (e.g., time_levs(1)state%h(:,:))
+ ! plus grid meta-data
+ ! Output: domain - upon exit, time level 2 (e.g., time_levs(2)%state%h(:,:)) contains
+ ! model state advanced forward in time by dt seconds
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+ implicit none
+
+ type (domain_type), intent(inout) :: domain
+ real (kind=RKIND), intent(in) :: dt
+
+ type (block_type), pointer :: block
+ type (state_type), pointer :: state
+ 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 :: uReconstructX, uReconstructY, uReconstructZ, &
+ uReconstructZonal, uReconstructMeridional
+
+ integer, dimension(:), allocatable :: mask
+
+ integer :: nVertLevels, k
+
+ integer :: err
+
+
+ block => domain % blocklist
+ do while (associated(block))
+ 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
+
+ uReconstructX => state % uReconstructX % array
+ uReconstructY => state % uReconstructY % array
+ uReconstructZ => state % uReconstructZ % array
+ uReconstructZonal => state % uReconstructZonal % array
+ uReconstructMeridional => state % uReconstructMeridional % array
+
+ !layerThickness_tend => block % tend % layerThickness % array
+ thickness_tend => block % tend % thickness % array
+
+ k = size(thickness)
+ allocate(mask(k))
+
+
+ ! 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
+ !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
+
+ !Optionally print some information about the new state
+ print *, 'thickness maxval:', maxval(thickness)
+ mask = 0
+ where (thickness .lt. 0.0)
+ mask = 1
+ thickness = 0.0
+ end where
+ print *,'Cells with negative thickness (set to 0):',sum(mask)
+
+ mask = 0
+ where (thickness .gt. 0.0)
+ mask = 1
+ end where
+ print *,'Cells with nonzero thickness:', sum(mask)
+
+ ! call land_ice_diagnostic_solve()
+
+ ! Reconstruct velocities for this time step
+ call mpas_reconstruct(mesh, normalVelocity, &
+ uReconstructX, uReconstructY, uReconstructZ, &
+ uReconstructZonal, uReconstructMeridional)
+
+ block => block % next
+ end do
+
+ deallocate(mask)
+
+ end subroutine land_ice_time_integrator_ForwardEuler
+
+
+
+end module land_ice_time_integration_ForwardEuler
+
+
Modified: branches/land_ice_projects/implement_core/src/framework/mpas_constants.F
===================================================================
--- branches/land_ice_projects/implement_core/src/framework/mpas_constants.F 2012-04-17 18:59:11 UTC (rev 1789)
+++ branches/land_ice_projects/implement_core/src/framework/mpas_constants.F 2012-04-17 20:15:32 UTC (rev 1790)
@@ -11,6 +11,7 @@
real (kind=RKIND), parameter :: cv = 716. ! cp - rgas
real (kind=RKIND), parameter :: cvpm = -.71385842 ! -cv/cp
real (kind=RKIND), parameter :: prandtl = 1.0
+ real (kind=RKIND), parameter :: SecondsInYear = 31536000.0 ! For 365 day year
contains
</font>
</pre>