<p><b>mhoffman@lanl.gov</b> 2012-03-20 11:13:33 -0600 (Tue, 20 Mar 2012)</p><p>BRANCH COMMIT<br>
Created a Shallow Ice Approximation velocity solver within MPAS. <br>
This will be useful for testing advection and other components without needing to compile with (and solve with) LifeV. The velocity solver can be chosen in the namelist file. I tested it and it matches well with velocities from CISM's SIA solver.<br>
</p><hr noshade><pre><font color="gray">Modified: branches/land_ice_projects/implement_core/src/core_land_ice/Makefile
===================================================================
--- branches/land_ice_projects/implement_core/src/core_land_ice/Makefile 2012-03-20 17:07:23 UTC (rev 1684)
+++ branches/land_ice_projects/implement_core/src/core_land_ice/Makefile 2012-03-20 17:13:33 UTC (rev 1685)
@@ -4,6 +4,7 @@
mpas_land_ice_time_integration.o \
mpas_land_ice_vel.o \
mpas_land_ice_lifev.o \
+ mpas_land_ice_sia.o \
mpas_land_ice_global_diagnostics.o
all: core_land_ice
@@ -13,10 +14,12 @@
mpas_land_ice_time_integration.o: mpas_land_ice_vel.o
-mpas_land_ice_vel.o: mpas_land_ice_lifev.o
+mpas_land_ice_vel.o: mpas_land_ice_lifev.o mpas_land_ice_sia.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
Modified: branches/land_ice_projects/implement_core/src/core_land_ice/Registry
===================================================================
--- branches/land_ice_projects/implement_core/src/core_land_ice/Registry 2012-03-20 17:07:23 UTC (rev 1684)
+++ branches/land_ice_projects/implement_core/src/core_land_ice/Registry 2012-03-20 17:13:33 UTC (rev 1685)
@@ -16,6 +16,8 @@
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 character io config_input_name grid.nc
namelist character io config_output_name output.nc
namelist character io config_restart_name restart.nc
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 2012-03-20 17:07:23 UTC (rev 1684)
+++ branches/land_ice_projects/implement_core/src/core_land_ice/mpas_land_ice_global_diagnostics.F 2012-03-20 17:13:33 UTC (rev 1685)
@@ -14,7 +14,7 @@
subroutine land_ice_compute_global_diagnostics(dminfo, state, grid, timeIndex, dt)
! Note: this routine assumes that there is only one block per processor. No looping
- ! is preformed over blocks.
+ ! is performed over blocks.
! dminfo is the domain info needed for global communication
! state contains the state variables needed to compute global diagnostics
! grid conains the meta data about the grid
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-03-20 17:07:23 UTC (rev 1684)
+++ branches/land_ice_projects/implement_core/src/core_land_ice/mpas_land_ice_lifev.F 2012-03-20 17:13:33 UTC (rev 1685)
@@ -314,8 +314,6 @@
integer :: index_temperature
-
-
err = 0
thickness => state % thickness % array
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-03-20 17:07:23 UTC (rev 1684)
+++ branches/land_ice_projects/implement_core/src/core_land_ice/mpas_land_ice_mpas_core.F 2012-03-20 17:13:33 UTC (rev 1685)
@@ -28,7 +28,7 @@
real (kind=RKIND) :: dt
type (block_type), pointer :: block
- integer :: err
+ integer :: i, err
err = 0
@@ -46,11 +46,21 @@
do while (associated(block))
call mpas_init_block(block, block % mesh, dt)
block % state % time_levs(1) % state % xtime % scalar = startTimeStamp
+ ! Make sure all time levels have a copy of the initial state prior the first velocity solve.
+ do i=2,nTimeLevs
+ call mpas_copy_state(block % state % time_levs(i) % state, block % state % time_levs(1) % state)
+ end do
+
block => block % next
end do
current_outfile_frames = 0
+ if(err == 1) then
+ print *, "An error has occurred in mpas_core_init. Aborting..."
+ call mpas_dmpar_global_abort("An error has occurred in mpas_core_init. Aborting...")
+ endif
+
end subroutine mpas_core_init
@@ -76,7 +86,7 @@
if (trim(config_stop_time) /= "none") then
call mpas_set_time(curr_time=stopTime, dateTimeString=config_stop_time, ierr=ierr)
if(startTime + runduration /= stopTime) then
- write(0,*) 'Warning: config_run_duration and config_stop_time are inconsitent: using config_run_duration.'
+ write(0,*) 'Warning: config_run_duration and config_stop_time are inconsistent: using config_run_duration.'
end if
end if
else if (trim(config_stop_time) /= "none") then
Added: 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 (rev 0)
+++ branches/land_ice_projects/implement_core/src/core_land_ice/mpas_land_ice_sia.F 2012-03-20 17:13:33 UTC (rev 1685)
@@ -0,0 +1,438 @@
+!|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
+!
+! land_ice_sia
+!
+!> \MPAS land-ice SIA velocity driver
+!> \author Matt Hoffman
+!> \date 16 March 2012
+!> \version SVN:$Id:$
+!> \details
+!> This module contains the routines for calculating velocity using the shallow ice approximation.
+!>
+!
+!-----------------------------------------------------------------------
+
+module land_ice_sia
+
+ use mpas_grid_types
+ use mpas_configure
+ use mpas_dmpar
+
+ implicit none
+ private
+ save
+
+ !--------------------------------------------------------------------
+ !
+ ! Public parameters
+ !
+ !--------------------------------------------------------------------
+
+ !--------------------------------------------------------------------
+ !
+ ! Public member functions
+ !
+ !--------------------------------------------------------------------
+ public :: land_ice_sia_init, &
+ land_ice_sia_finalize, &
+ land_ice_sia_block_init, &
+ land_ice_sia_solve
+
+ !--------------------------------------------------------------------
+ !
+ ! Private module variables
+ !
+ !--------------------------------------------------------------------
+
+ ! halo exchange arrays
+ integer, dimension(:), pointer :: sendCellsArray, &
+ recvCellsArray, &
+ sendVerticesArray, &
+ recvVerticesArray, &
+ sendEdgesArray, &
+ recvEdgesArray
+
+!***********************************************************************
+
+contains
+
+!***********************************************************************
+!
+! routine land_ice_sia_init
+!
+!> \brief Initializes velocity solver
+!> \author Matt Hoffman/Xylar Asay-Davis
+!> \date 16 March 2012
+!> \version SVN:$Id$
+!> \details
+!> This routine initializes the ice velocity solver.
+!
+!-----------------------------------------------------------------------
+
+ subroutine land_ice_sia_init(domain, err)
+
+ !-----------------------------------------------------------------
+ !
+ ! input variables
+ !
+ !-----------------------------------------------------------------
+
+ type (domain_type), intent(inout) :: domain
+
+ !-----------------------------------------------------------------
+ !
+ ! input/output variables
+ !
+ !-----------------------------------------------------------------
+
+ !-----------------------------------------------------------------
+ !
+ ! output variables
+ !
+ !-----------------------------------------------------------------
+
+ integer, intent(out) :: err !< Output: error flag
+
+ !-----------------------------------------------------------------
+ !
+ ! local variables
+ !
+ !-----------------------------------------------------------------
+
+ integer :: nCells, nEdges, nVertices, nCellsSolve, nEdgesSolve, nVerticesSolve, nVertLevels
+
+ integer, dimension(:,:), pointer :: cellsOnEdge, cellsOnVertex, verticesOnCell, verticesOnEdge
+
+ real (kind=RKIND), dimension(:), pointer :: xCell, yCell, zCell
+
+ logical, dimension(:), pointer :: verticesMask
+
+ logical :: keepVertex
+
+ err = 0
+
+ ! Note: the following code assumes that blocklist contains a single block
+ ! Significant modification will be needed when multiple blocks per processor
+ ! are supported.
+
+ !extract data from domain
+ nCellsSolve = domain % blocklist % mesh % nCellsSolve
+ nEdgesSolve = domain % blocklist % mesh % nEdgesSolve
+ nVerticesSolve = domain % blocklist % mesh % nVerticesSolve
+
+ nCells = domain % blocklist % mesh % nCells
+ nEdges = domain % blocklist % mesh % nEdges
+ nVertices = domain % blocklist % mesh % nVertices
+ nVertLevels = domain % blocklist % mesh % nVertLevels
+
+ cellsOnEdge => domain % blocklist % mesh % cellsOnEdge % array
+ cellsOnVertex => domain % blocklist % mesh % cellsOnVertex % array
+ verticesOnCell => domain % blocklist % mesh % verticesOnCell % array
+ verticesOnEdge => domain % blocklist % mesh % verticesOnEdge % array
+
+ xCell => domain % blocklist % mesh % xCell % array
+ yCell => domain % blocklist % mesh % yCell % array
+ zCell => domain % blocklist % mesh % zCell % array
+
+ !allocate(verticesMask(nVertices))
+
+ !verticesMask = .true. ! all vertices are part of the domain, since thickness isn't necessarily valid here
+
+
+ !build send and receive arrays using exchange_list
+
+ call array_from_exchange_list(domain % blocklist % parinfo % verticesToSend, sendVerticesArray)
+ call array_from_exchange_list(domain % blocklist % parinfo % verticesToRecv, recvVerticesArray)
+ call array_from_exchange_list(domain % blocklist % parinfo % cellsToSend, sendCellsArray)
+ call array_from_exchange_list(domain % blocklist % parinfo % cellsToRecv, recvCellsArray)
+ call array_from_exchange_list(domain % blocklist % parinfo % edgesToSend, sendEdgesArray)
+ call array_from_exchange_list(domain % blocklist % parinfo % edgesToRecv, recvEdgesArray)
+
+
+ !--------------------------------------------------------------------
+
+ end subroutine land_ice_sia_init
+
+!***********************************************************************
+!
+! routine land_ice_sia_finalize
+!
+!> \brief finalizes velocity solver
+!> \author Matt Hoffman/Xylar Asay-Davis
+!> \date 16 March 2012
+!> \version SVN:$Id$
+!> \details
+!> This routine initializes the ice velocity solver.
+!
+!-----------------------------------------------------------------------
+
+ subroutine land_ice_sia_finalize(domain, err)
+
+ !-----------------------------------------------------------------
+ !
+ ! input variables
+ !
+ !-----------------------------------------------------------------
+
+ type (domain_type), intent(inout) :: domain
+
+ !-----------------------------------------------------------------
+ !
+ ! input/output variables
+ !
+ !-----------------------------------------------------------------
+
+ !-----------------------------------------------------------------
+ !
+ ! output variables
+ !
+ !-----------------------------------------------------------------
+
+ integer, intent(out) :: err !< Output: error flag
+
+ !-----------------------------------------------------------------
+ !
+ ! local variables
+ !
+ !-----------------------------------------------------------------
+
+ err = 0
+
+ deallocate(sendCellsArray, &
+ recvCellsArray, &
+ sendVerticesArray, &
+ recvVerticesArray, &
+ sendEdgesArray, &
+ recvEdgesArray)
+
+ !--------------------------------------------------------------------
+
+ end subroutine land_ice_sia_finalize
+
+!***********************************************************************
+!
+! routine land_ice_sia_block_init
+!
+!> \brief Initializes velocity solver
+!> \author Matt Hoffman/Xylar Asay-Davis
+!> \date 16 March 2012
+!> \version SVN:$Id$
+!> \details
+!> This routine initializes each block of the ice velocity solver.
+!
+!-----------------------------------------------------------------------
+
+ subroutine land_ice_sia_block_init(mesh, err)
+
+ !-----------------------------------------------------------------
+ !
+ ! input variables
+ !
+ !-----------------------------------------------------------------
+
+ type (mesh_type), intent(in) :: &
+ mesh !< Input: mesh information
+
+ !-----------------------------------------------------------------
+ !
+ ! input/output variables
+ !
+ !-----------------------------------------------------------------
+
+ !-----------------------------------------------------------------
+ !
+ ! output variables
+ !
+ !-----------------------------------------------------------------
+
+ integer, intent(out) :: err !< Output: error flag
+
+ !-----------------------------------------------------------------
+ !
+ ! local variables
+ !
+ !-----------------------------------------------------------------
+
+ err = 0
+
+ !--------------------------------------------------------------------
+
+ end subroutine land_ice_sia_block_init
+
+!***********************************************************************
+
+ subroutine land_ice_sia_solve(mesh, state, err)
+ use mpas_constants, only: gravity
+
+ !-----------------------------------------------------------------
+ !
+ ! input variables
+ !
+ !-----------------------------------------------------------------
+
+ type (mesh_type), intent(in) :: &
+ mesh !< Input: mesh information
+
+ !-----------------------------------------------------------------
+ !
+ ! input/output variables
+ !
+ !-----------------------------------------------------------------
+
+ type (state_type), intent(in) :: &
+ state !< Input/Output: state information
+
+ !-----------------------------------------------------------------
+ !
+ ! output variables
+ !
+ !-----------------------------------------------------------------
+
+ integer, intent(out) :: err !< Output: error flag
+
+ !-----------------------------------------------------------------
+ !
+ ! local variables
+ !
+ !-----------------------------------------------------------------
+
+ real (kind=RKIND), dimension(:), pointer :: thickness, layerThicknessFractions
+ 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
+ logical edgeIsIce
+ real rhoi !!!! THIS SHOULD BE A PHYSICAL PARAMETER ELSEWHERE
+
+ err = 0
+
+ ! Set needed variables and pointers
+ nCells = mesh % nCells
+ nEdges = mesh % nEdges
+ nVertLevels = mesh % nVertLevels
+
+ dcEdge => mesh % dcEdge % array
+ cellsOnEdge => mesh % cellsOnEdge % array
+ layerThicknessFractions => mesh % layerThicknessFractions % array
+
+ thickness => state % thickness % array
+ layerThickness => state % layerThickness % array
+ normalVelocity => state % normalVelocity % array
+
+ rhoi = config_ice_density
+
+ basalVelocity = 0.0 ! Assume no sliding
+
+ ratefactor = 1.0e-16
+
+ ! Determine the height of each layer above the bed
+ layerFractionAboveBed(1) = 0.5 * layerThicknessFractions(1)
+ do iLevel = 2, nVertLevels
+ layerFractionAboveBed(iLevel) = layerFractionAboveBed(iLevel - 1) + &
+ 0.5 * layerThicknessFractions(iLevel - 1) + 0.5 * layerThicknessFractions(iLevel)
+ end do
+
+ ! Loop over edges
+ do iEdge = 1, nEdges
+ ! Determine if ice
+ call is_edge_ice(thickness, cellsOnEdge, iEdge, edgeIsIce)
+ ! Only calculate velocity for edges that are part of the ice sheet.
+ ! This should be a mask calculated external this module.
+ ! Also, the velocity calculation should be valid for non-ice edges (i.e. returns 0).
+ if (edgeIsIce .eqv. .true.) then
+ ! Calculate slope, thickness at edge
+ ! These could/should be calculated externally to this subroutine
+ slopeOnEdge = (thickness(cellsOnEdge(1,iEdge)) - thickness(cellsOnEdge(2,iEdge)) ) &
+ / dcEdge(iEdge)
+ thicknessOnEdge = 0.5 * (thickness(cellsOnEdge(1,iEdge)) + thickness(cellsOnEdge(2,iEdge)))
+ ! Loop over layers
+ do iLevel = 1, nVertLevels
+ ! Calculate ratefactor (A) at edge - This should be calculated here as a function of temperature
+ ! ratefactor = f(T)
+ ! Calculate SIA velocity
+ normalVelocity(iLevel,iEdge) = basalVelocity + &
+ 0.5 * ratefactor * (rhoi * gravity)**3 * slopeOnEdge**3 * &
+ (thicknessOnEdge**4 - (thicknessOnEdge - thicknessOnEdge * layerFractionAboveBed(iLevel))**4)
+ end do ! Levels
+ endif
+ end do ! edges
+
+ !--------------------------------------------------------------------
+
+ end subroutine land_ice_sia_solve
+
+
+ ! private subroutines
+
+ subroutine is_edge_ice(thickness, cellsOnEdge, iEdge, edgeIsIce)
+ ! Determine if an edge is part of the ice sheet.
+ ! This should probably be calculated external to this module as a mask.
+ ! (and could be more sophisticated than this)
+
+ implicit none
+
+ real (kind=RKIND), dimension(:), intent(in) :: thickness
+ integer, dimension(:,:), intent(in) :: cellsOnEdge
+ integer, intent(in) :: iEdge
+ logical, intent(out) :: edgeIsIce
+
+ ! Local variables
+ integer :: cell1, cell2
+
+ cell1 = cellsOnEdge(1,iEdge)
+ cell2 = cellsOnEdge(2,iEdge)
+ ! Check if an edge has at least one adjacent cell with ice
+ if ( (thickness(cell1).gt.0.0) .or. (thickness(cell2).gt.0.0) ) then
+ edgeIsIce = .true.
+ else
+ edgeIsIce = .false.
+ endif
+
+ end subroutine is_edge_ice
+
+ !--------------------------------------------------------------------
+
+ subroutine array_from_exchange_list(exlist, array)
+ ! This is identical to the version in land_ice_lifev and should probably put in a common location.
+
+ implicit none
+
+ type (exchange_list), pointer, intent(in) :: exlist
+ type (exchange_list), pointer :: listPtr
+
+ integer, dimension(:), pointer, intent(out) :: array
+ integer :: dataSize, offset, i
+
+ dataSize = 1 !in first position we will store the size of the vector
+ listPtr => exlist
+ do while (associated(listPtr))
+ dataSize = dataSize + (listPtr % nlist) + 2
+ listPtr => listPtr % next
+ end do
+
+ allocate(array(dataSize))
+
+ array(1) = dataSize;
+ offset = 2;
+ listPtr => exlist
+ do while (associated(listPtr))
+ array(offset) = listPtr % procID
+ offset = offset + 1
+ array(offset) = listPtr % nlist
+ do i=1,listPtr % nlist
+ array(i+offset) = listPtr % list(i) -1
+ end do
+ offset = offset + listPtr % nlist + 1
+ listPtr => listPtr % next
+ end do
+
+ end subroutine array_from_exchange_list
+
+!***********************************************************************
+
+end module land_ice_sia
+
+!|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
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-03-20 17:07:23 UTC (rev 1684)
+++ branches/land_ice_projects/implement_core/src/core_land_ice/mpas_land_ice_time_integration.F 2012-03-20 17:13:33 UTC (rev 1685)
@@ -26,11 +26,13 @@
character(len=*), intent(in) :: timeStamp
type (block_type), pointer :: block
+ type (state_type), pointer :: state
- real (kind=RKIND), dimension(:), pointer :: thickness, layerThicknessFractions
- real (kind=RKIND), dimension(:,:), pointer :: layerThickness
+ real (kind=RKIND), dimension(:,:), pointer :: normalVelocity, uReconstructX, &
+ uReconstructY, uReconstructZ, uReconstructZonal, uReconstructMeridional
- integer :: nVertLevels, nCells, iCell, iLevel, err
+ integer :: err
+
! if (trim(config_time_integration) == 'RK4') then
! call land_ice_rk4(domain, dt)
@@ -43,26 +45,23 @@
! a temporary test that calls the velocity solver on the state at the current time
block => domain % blocklist
- nCells = block % mesh % nCells
- nVertLevels = block % mesh % nVertLevels
- layerThicknessFractions => block % mesh % layerThicknessFractions % array
+ 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
- thickness => block % state % time_levs(1) % state % thickness % array
- layerThickness => block % state % time_levs(1) % state % layerThickness % array
+ call land_ice_vel_solve(block % mesh, state, err)
+ state % xtime % scalar = timeStamp
- ! set up layerThickness from thickness using the layerThicknessFractions
- do iCell = 1, nCells
- do iLevel = 1, nVertLevels
- layerThickness(iLevel,iCell) = layerThicknessFractions(iLevel)*thickness(iCell)
- end do
- end do
+ ! 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 => domain % blocklist
- call land_ice_vel_solve(block % mesh, block % state % time_levs(1) % state, err)
-
- block => domain % blocklist
- do while (associated(block))
- block % state % time_levs(2) % state % xtime % scalar = timeStamp
block => block % next
end do
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-03-20 17:07:23 UTC (rev 1684)
+++ branches/land_ice_projects/implement_core/src/core_land_ice/mpas_land_ice_vel.F 2012-03-20 17:13:33 UTC (rev 1685)
@@ -17,6 +17,7 @@
use mpas_grid_types
use mpas_configure
use land_ice_lifev
+ use land_ice_sia
implicit none
private
@@ -94,10 +95,20 @@
!-----------------------------------------------------------------
err = 0
+
+ write(*,*) 'Using ', trim(config_dycore), ' dynamical core.'
+ select case (config_dycore)
+ case ('L1L2')
+ call land_ice_lifev_init(domain, err)
+ case ('SIA')
+ call land_ice_sia_init(domain, err)
+ case default
+ write(*,*) trim(config_dycore), ' is not a valid land ice dycore option.'
+ err = 1
+ return
+ end select
- call land_ice_lifev_init(domain, err)
-
!--------------------------------------------------------------------
end subroutine land_ice_vel_init
@@ -147,8 +158,18 @@
err = 0
- call land_ice_lifev_finalize(domain, err)
+ select case (config_dycore)
+ case ('L1L2')
+ call land_ice_lifev_finalize(domain, err)
+ case ('SIA')
+ call land_ice_sia_finalize(domain, err)
+ case default
+ write(*,*) trim(config_dycore), ' is not a valid land ice dycore option.'
+ err = 1
+ return
+ end select
+
!--------------------------------------------------------------------
end subroutine land_ice_vel_finalize
@@ -199,8 +220,18 @@
err = 0
- call land_ice_lifev_block_init(mesh, err)
+ select case (config_dycore)
+ case ('L1L2')
+ call land_ice_lifev_block_init(mesh, err)
+ case ('SIA')
+ call land_ice_sia_block_init(mesh, err)
+ case default
+ write(*,*) trim(config_dycore), ' is not a valid land ice dycore option.'
+ err = 1
+ return
+ end select
+
!--------------------------------------------------------------------
end subroutine land_ice_vel_block_init
@@ -209,6 +240,8 @@
subroutine land_ice_vel_solve(mesh, state, err)
+ use land_ice_SIA
+
!-----------------------------------------------------------------
!
! input variables
@@ -243,7 +276,16 @@
err = 0
- call land_ice_lifev_solve(mesh, state, err)
+ select case (config_dycore)
+ case ('L1L2')
+ call land_ice_lifev_solve(mesh, state, err)
+ case ('SIA')
+ call land_ice_sia_solve(mesh, state, err)
+ case default
+ write(*,*) trim(config_dycore), ' is not a valid land ice dycore option.'
+ err = 1
+ return
+ end select
!--------------------------------------------------------------------
</font>
</pre>