<p><b>mhoffman@lanl.gov</b> 2012-04-26 12:45:50 -0600 (Thu, 26 Apr 2012)</p><p>BRANCH COMMIT<br>
<br>
Added logic for, and calls to, FO and Stokes velocity solvers in LifeV.<br>
Also, general cleanup and reorganization of the land_ice_lifev module (e.g. added the vertexMask to Registry and moved vertexMask calculation to the land_ice_tendency module, moved most of the init to block_init).<br>
I've tested the changes with the Fortran SIA velocity solver.<br>
I've tested the LifeV implementation as much as I can without having working LifeV libraries (i.e. the land_ice_lifev module compiles LIFEV=true, land_ice_model.exe runs successfully with LIFEV=false but with config_dycore='L1L2').<br>
</p><hr noshade><pre><font color="gray">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-26 18:44:42 UTC (rev 1818)
+++ branches/land_ice_projects/implement_core/src/core_land_ice/Registry 2012-04-26 18:45:50 UTC (rev 1819)
@@ -12,7 +12,7 @@
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.
+% valid dycores are SIA, L1L2, FO, 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
@@ -57,6 +57,9 @@
%
var persistent text xtime ( Time ) 2 ro xtime state - -
+
+% MESH VARIABLES =====================
+
var persistent real latCell ( nCells ) 0 iro latCell mesh - -
var persistent real lonCell ( nCells ) 0 iro lonCell mesh - -
var persistent real xCell ( nCells ) 0 iro xCell mesh - -
@@ -123,6 +126,9 @@
var persistent real layerCenterSigma ( nVertLevels ) 0 - layerCenterSigma mesh - -
var persistent real layerInterfaceSigma ( nVertLevelsPlus1 ) 0 - layerInterfaceSigma mesh - -
+
+% STATE VARIABLES =====================
+
% Prognostic variables: read from input, saved in restart, and written to output
var persistent real thickness ( nCells Time ) 2 iro thickness state - -
var persistent real bedTopography ( nCells Time ) 2 iro bedTopography state - -
@@ -152,5 +158,13 @@
var persistent real layerThicknessVertex ( nVertLevels nVertices Time ) 2 - layerThicknessVertex state - -
var persistent real layerThicknessEdge ( nVertLevels nEdges Time ) 2 - layerThicknessEdge state - -
+% Masks: calculated diagnostically at each time step
+var persistent integer cellMask ( nCells Time ) 2 o cellMask state - -
+var persistent integer edgeMask ( nEdges Time ) 2 o edgeMask state - -
+var persistent integer vertexMask ( nVertices Time ) 2 o vertexMask state - -
+% Needed for LifeV solver to know if the region to solve on the block's domain has changed (treated as a logical):
+% THIS SHOULD BE CHANGED TO INTEGER TYPE ONCE THE TRUNK IS MERGED TO THIS BRANCH AND field0dinteger IS DEFINED IN framework/mpas_grid_types.F
+var persistent real vertexMaskChanged ( ) 2 - vertexMaskChanged state - -
+
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-26 18:44:42 UTC (rev 1818)
+++ branches/land_ice_projects/implement_core/src/core_land_ice/mpas_land_ice_lifev.F 2012-04-26 18:45:50 UTC (rev 1819)
@@ -40,7 +40,13 @@
land_ice_lifev_finalize, &
land_ice_lifev_block_init, &
land_ice_lifev_solve
+ !--------------------------------------------------------------------
+ !
+ ! Private module variables
+ !
+ !--------------------------------------------------------------------
+
!***********************************************************************
contains
@@ -88,83 +94,16 @@
!
!-----------------------------------------------------------------
- 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
-
- ! halo exchange arrays
- integer, dimension(:), pointer :: sendCellsArray, &
- recvCellsArray, &
- sendVerticesArray, &
- recvVerticesArray, &
- sendEdgesArray, &
- recvEdgesArray
-
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.
+#ifdef USE_LIFEV
+ ! These calls are needed for using any of the LifeV velocity solvers
- !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
-
- !call lifev_init(dminfo % comm, nVerteces)
-
- !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)
-
-#ifdef USE_LIFEV
!call lifeV and set the grid of the velocity solver
call velocity_solver_init_mpi(domain % dminfo % comm)
- !zCell is supposed to be zero for L1L2 and FO solvers
- !nVertLevels should be equal to nVertLevelsSolve (no split the domain in the vertical direction)
- call velocity_solver_set_grid_data(nCells, nEdges, nVertices, nVertLevels, nCellsSolve, nEdgesSolve, nVerticesSolve, &
- cellsOnEdge, cellsOnVertex, verticesOnCell, verticesOnEdge, xCell, yCell, zCell, &
- sendCellsArray, recvCellsArray, sendEdgesArray, recvEdgesArray, sendVerticesArray, recvVerticesArray)
#endif
- !these ca be deallocated because they have been copied on the c++ side
- deallocate(sendCellsArray, &
- recvCellsArray, &
- sendVerticesArray, &
- recvVerticesArray, &
- sendEdgesArray, &
- recvEdgesArray)
-
!--------------------------------------------------------------------
end subroutine land_ice_lifev_init
@@ -215,6 +154,7 @@
err = 0
#ifdef USE_LIFEV
+ ! This call is needed for using any of the LifeV velocity solvers
call velocity_solver_finalize()
#endif
@@ -235,7 +175,7 @@
!
!-----------------------------------------------------------------------
- subroutine land_ice_lifev_block_init(mesh, err)
+ subroutine land_ice_lifev_block_init(block, err)
!-----------------------------------------------------------------
!
@@ -243,8 +183,8 @@
!
!-----------------------------------------------------------------
- type (mesh_type), intent(in) :: &
- mesh !< Input: mesh information
+ type (block_type), intent(in) :: &
+ block !< Input: mesh information
!-----------------------------------------------------------------
!
@@ -266,8 +206,73 @@
!
!-----------------------------------------------------------------
+ integer :: nCells, nEdges, nVertices, nCellsSolve, nEdgesSolve, nVerticesSolve, nVertLevels
+
+ integer, dimension(:,:), pointer :: cellsOnEdge, cellsOnVertex, verticesOnCell, verticesOnEdge
+
+ real (kind=RKIND), dimension(:), pointer :: xCell, yCell, zCell
+
+ ! halo exchange arrays
+ integer, dimension(:), pointer :: sendCellsArray, &
+ recvCellsArray, &
+ sendVerticesArray, &
+ recvVerticesArray, &
+ sendEdgesArray, &
+ recvEdgesArray
+
err = 0
+ ! Note: LifeV only supports one block per processor, but this has (hopefully)
+ ! been written to work if that were to change. (That's why all these LifeV init
+ ! calls are in land_ice_lifev_block_init instead of land_ice_lifev_init.)
+
+ !extract data from domain
+ nCellsSolve = block % mesh % nCellsSolve
+ nEdgesSolve = block % mesh % nEdgesSolve
+ nVerticesSolve = block % mesh % nVerticesSolve
+
+ nCells = block % mesh % nCells
+ nEdges = block % mesh % nEdges
+ nVertices = block % mesh % nVertices
+ nVertLevels = block % mesh % nVertLevels
+
+ cellsOnEdge => block % mesh % cellsOnEdge % array
+ cellsOnVertex => block % mesh % cellsOnVertex % array
+ verticesOnCell => block % mesh % verticesOnCell % array
+ verticesOnEdge => block % mesh % verticesOnEdge % array
+
+ xCell => block % mesh % xCell % array
+ yCell => block % mesh % yCell % array
+ zCell => block % mesh % zCell % array
+
+
+ !build send and receive arrays using exchange_list
+ call array_from_exchange_list(block % parinfo % verticesToSend, sendVerticesArray)
+ call array_from_exchange_list(block % parinfo % verticesToRecv, recvVerticesArray)
+ call array_from_exchange_list(block % parinfo % cellsToSend, sendCellsArray)
+ call array_from_exchange_list(block % parinfo % cellsToRecv, recvCellsArray)
+ call array_from_exchange_list(block % parinfo % edgesToSend, sendEdgesArray)
+ call array_from_exchange_list(block % parinfo % edgesToRecv, recvEdgesArray)
+
+#ifdef USE_LIFEV
+ ! These calls are needed for using any of the LifeV velocity solvers
+
+ !zCell is supposed to be zero for L1L2 and FO solvers
+ !nVertLevels should be equal to nVertLevelsSolve (no split the domain in the vertical direction)
+ call velocity_solver_set_grid_data(nCells, nEdges, nVertices, nVertLevels, nCellsSolve, nEdgesSolve, nVerticesSolve, &
+ cellsOnEdge, cellsOnVertex, verticesOnCell, verticesOnEdge, xCell, yCell, zCell, &
+ sendCellsArray, recvCellsArray, sendEdgesArray, recvEdgesArray, sendVerticesArray, recvVerticesArray)
+#endif
+
+ !these can be deallocated because they have been copied on the c++ side
+ deallocate(sendCellsArray, &
+ recvCellsArray, &
+ sendVerticesArray, &
+ recvVerticesArray, &
+ sendEdgesArray, &
+ recvEdgesArray)
+
+
!--------------------------------------------------------------------
end subroutine land_ice_lifev_block_init
@@ -308,74 +313,66 @@
!
!-----------------------------------------------------------------
real (kind=RKIND), dimension(:), pointer :: &
- thickness, lowerSurface, upperSurface, beta, levelsRatio
+ thickness, lowerSurface, upperSurface, beta, LayerThicknessFractions
real (kind=RKIND), dimension(:,:), pointer :: &
normalVelocity
real (kind=RKIND), dimension(:,:,:), pointer :: &
tracers
+ integer, dimension(:), pointer :: vertexMask
- !temporary variables
- logical, dimension(:), pointer :: verticesMask
- logical :: keep_vertex
- integer :: nVertices, i, k, iCell, nVertLevels
-
-
-
-
integer :: index_temperature
+ real (kind=RKIND) :: vertexMaskChanged ! THIS TYPE NEEDS TO BE CHANGE TO INTEGER ONCE REGISTRY IS UPDATED
err = 0
+
+ ! Mesh variables
+ LayerThicknessFractions => mesh % LayerThicknessFractions % array
+ ! State variables
normalVelocity => state % normalVelocity % array
thickness => state % thickness % array
lowerSurface => state % lowerSurface % array
upperSurface => state % upperSurface % array
beta => state % beta % array
tracers => state % tracers % array
-
index_temperature = state % index_temperature
+ vertexMask => state % vertexMask % array
+ vertexMaskChanged = state % vertexMaskChanged % scalar
- nVertices = mesh % nVertices
- allocate(verticesMask(nVertices))
- do i = 1,nVertices
- keep_vertex = .false.
- do k = 1, 3
- iCell = mesh % cellsOnVertex % array(k,i)
- keep_vertex = keep_vertex.or.(thickness(iCell) > 0.01)
- end do
- verticesMask(i) = keep_vertex
- end do
-
- nVertLevels = mesh % nVertLevels
- allocate(levelsRatio(nVertLevels))
- do i = 1,nVertLevels
- levelsRatio(i) = 1./real(nVertLevels);
- end do
-
#ifdef USE_LIFEV
- !<---- calls to be made only when mask changes
- call velocity_solver_compute_2d_grid(verticesMask)
+ ! LifeV calls to be made only when vertex mask changes
+ if (vertexMaskChanged .eq. 1) then
+ call velocity_solver_compute_2d_grid(vertexMask)
+ select case (config_dycore)
+ case ('L1L2')
+ call velocity_solver_init_L1L2(LayerThicknessFractions)
+ case ('FO')
+ call velocity_solver_extrude_3d_grid(LayerThicknessFractions, lowerSurface, thickness)
+ call velocity_solver_init_FO(LayerThicknessFractions)
+ case ('Stokes')
+ call velocity_solver_extrude_3d_grid(LayerThicknessFractions, lowerSurface, thickness)
+ call valocity_solver_init_stokes(LayerThicknessFractions)
+ end select
+ endif
- call velocity_solver_init_L1L2(levelsRatio)
-
- !if FO or Stokes we need this
- !call velocity_solver_extrude_3d_grid(levelsRatio, lowerSurface, thickness)
- !call velocity_solver_init_FO(levelsRatio)
- !call valocity_solver_init_stokes(levelsRatio)
+ ! LifeV calls to be made every time step (solve velocity!)
+ select case (config_dycore)
+ case ('L1L2')
+ call velocity_solver_solve_L1L2(lowerSurface, thickness, beta, tracers(index_temperature,:,:), normalVelocity, normalVelocity)
+ ! Optional calls to have LifeV output data files
+ call velocity_solver_export_2d_data(lowerSurface, thickness, beta)
+ call velocity_solver_export_L1L2_velocity();
+ case ('FO')
+ call velocity_solver_solve_FO(lowerSurface, thickness, beta, tracers(index_temperature,:,:), normalVelocity, normalVelocity)
+ call velocity_solver_export_2d_data(lowerSurface, thickness, beta)
+ call velocity_solver_export_FO_velocity();
+ case ('Stokes')
+ ! call stokes
+ end select
- !-->
-
- call velocity_solver_solve_L1L2(lowerSurface, thickness, beta, tracers(index_temperature,:,:), normalVelocity, normalVelocity)
- !call velocity_solver_solve_FO(lowerSurface, thickness, beta, tracers(index_temperature,:,:), normalVelocity, normalVelocity)
- call velocity_solver_export_L1L2_velocity();
- call velocity_solver_export_2d_data(lowerSurface, thickness, beta)
- !call velocity_solver_solve_FO(lowerSurface, thckdata, betadata, temperaturedata, u, v)
- !call velocity_solver_export_FO_velocity();
-
-
#endif
!--------------------------------------------------------------------
@@ -383,6 +380,10 @@
end subroutine land_ice_lifev_solve
+!***********************************************************************
+!***********************************************************************
+
+
! private subroutines
subroutine array_from_exchange_list(exlist, 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-04-26 18:44:42 UTC (rev 1818)
+++ branches/land_ice_projects/implement_core/src/core_land_ice/mpas_land_ice_mpas_core.F 2012-04-26 18:45:50 UTC (rev 1819)
@@ -53,6 +53,7 @@
block => domain % blocklist
do while (associated(block))
+ ! Initialize blocks
call mpas_init_block(block, block % mesh, dt)
! Assign initial time stamp
block % state % time_levs(1) % state % xtime % scalar = startTimeStamp
@@ -164,9 +165,10 @@
call mpas_rbf_interp_initialize(mesh)
call mpas_init_reconstruct(mesh)
- call land_ice_vel_block_init(mesh, err_tmp)
+ call land_ice_vel_block_init(block, err_tmp)
err = ior(err, err_tmp)
+
! Initialize state ====
call lice_diagnostic_solve(mesh, state, err)
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-26 18:44:42 UTC (rev 1818)
+++ branches/land_ice_projects/implement_core/src/core_land_ice/mpas_land_ice_sia.F 2012-04-26 18:45:50 UTC (rev 1819)
@@ -222,7 +222,7 @@
!
!-----------------------------------------------------------------------
- subroutine land_ice_sia_block_init(mesh, err)
+ subroutine land_ice_sia_block_init(block, err)
!-----------------------------------------------------------------
!
@@ -230,8 +230,8 @@
!
!-----------------------------------------------------------------
- type (mesh_type), intent(in) :: &
- mesh !< Input: mesh information
+ type (block_type), intent(in) :: &
+ block !< Input: mesh information
!-----------------------------------------------------------------
!
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-04-26 18:44:42 UTC (rev 1818)
+++ branches/land_ice_projects/implement_core/src/core_land_ice/mpas_land_ice_tendency.F 2012-04-26 18:45:50 UTC (rev 1819)
@@ -313,21 +313,27 @@
! local variables
!
!-----------------------------------------------------------------
- integer :: iCell, iLevel, nCells, nVertLevels
+ integer :: iCell, iLevel, nCells, nVertLevels, nVertices
real (kind=RKIND), dimension(:), pointer :: thickness, upperSurface, &
lowerSurface, bedTopography, layerThicknessFractions
real (kind=RKIND), dimension(:,:), pointer :: layerThickness
+ integer, dimension(:), pointer :: vertexMask
+ integer, dimension(:,:), pointer :: cellsOnVertex
+ integer :: keep_vertex, i, k
nCells = mesh % nCells
+ nVertices = mesh % nVertices
nVertLevels = mesh % nVertLevels
layerThicknessFractions => mesh % layerThicknessFractions % array
+ cellsOnVertex => mesh % cellsOnVertex % array
thickness => state % thickness % array
upperSurface => state % upperSurface % array
lowerSurface => state % lowerSurface % array
bedTopography => state % bedTopography % array
layerThickness => state % layerThickness % array
+ vertexMask => state % vertexMask % array
! no ice shelves -- lower surface same as bed topography
@@ -347,8 +353,21 @@
!\todo Calculate h_edge that can be used by SIA solve and FO advection scheme. ocn_diagnostic_solve in mpas_ocn_tendency.F has 2, 3, & 4th order calculations for h_edge that can be used. Right now I am using 2nd order h_edge in SIA solve and 1st order h_edge in FO advection.
- !\todo Calculate masks. Right now sia driver and Forward Euler time integration have need of a mask and both do it on their own. We need to define how the masks should work.
+ !\todo Calculate masks. Right now sia driver and Forward Euler time integration have need of a mask and both do it on their own. We need to define how the masks should work and implement them here.
+ ! simple vertex mask needed by LifeV
+ do i = 1,nVertices
+ keep_vertex = 0
+ do k = 1, 3
+ iCell = cellsOnVertex(k,i)
+ if (thickness(iCell) .gt. 0.01) then
+ keep_vertex = 1
+ endif
+ end do
+ vertexMask(i) = keep_vertex
+ end do
+
+
end subroutine lice_diagnostic_solve
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-26 18:44:42 UTC (rev 1818)
+++ branches/land_ice_projects/implement_core/src/core_land_ice/mpas_land_ice_time_integration_forwardeuler.F 2012-04-26 18:45:50 UTC (rev 1819)
@@ -138,7 +138,14 @@
! === Calculate diagnostic variables for new state =====================
call lice_diagnostic_solve(mesh, stateNew, err) ! perhaps velocity solve should move in here.
+ ! Determine if the vertex mask changed during this time step (needed for LifeV)
+ if ( sum(stateNew % vertexMask % array - stateOld % vertexMask % array) .ne. 0 ) then
+ stateNew % vertexMaskChanged % scalar = 1
+ else
+ stateNew % vertexMaskChanged % scalar = 0
+ endif
+
! Compute the diagnostic velocity for this time step using the updated (new) state
! Pass in the old velocity to be used as an initial guess for velocity solvers that need one.
normalVelocityNew = normalVelocityOld
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-26 18:44:42 UTC (rev 1818)
+++ branches/land_ice_projects/implement_core/src/core_land_ice/mpas_land_ice_vel.F 2012-04-26 18:45:50 UTC (rev 1819)
@@ -98,10 +98,14 @@
write(*,*) 'Using ', trim(config_dycore), ' dynamical core.'
select case (config_dycore)
+ case ('SIA')
+ call land_ice_sia_init(domain, err)
case ('L1L2')
call land_ice_lifev_init(domain, err)
- case ('SIA')
- call land_ice_sia_init(domain, err)
+ case ('FO')
+ call land_ice_lifev_init(domain, err)
+ case ('Stokes')
+ call land_ice_lifev_init(domain, err)
case default
write(*,*) trim(config_dycore), ' is not a valid land ice dycore option.'
err = 1
@@ -159,10 +163,15 @@
err = 0
select case (config_dycore)
+ case ('SIA')
+ call land_ice_sia_finalize(domain, err)
case ('L1L2')
call land_ice_lifev_finalize(domain, err)
- case ('SIA')
- call land_ice_sia_finalize(domain, err)
+ case ('FO')
+ call land_ice_lifev_finalize(domain, err)
+ case ('Stokes')
+ call land_ice_lifev_finalize(domain, err)
+
case default
write(*,*) trim(config_dycore), ' is not a valid land ice dycore option.'
err = 1
@@ -187,7 +196,7 @@
!
!-----------------------------------------------------------------------
- subroutine land_ice_vel_block_init(mesh, err)
+ subroutine land_ice_vel_block_init(block, err)
!-----------------------------------------------------------------
!
@@ -195,8 +204,8 @@
!
!-----------------------------------------------------------------
- type (mesh_type), intent(in) :: &
- mesh !< Input: mesh information
+ type (block_type), intent(in) :: &
+ block !< Input: mesh information
!-----------------------------------------------------------------
!
@@ -221,10 +230,14 @@
err = 0
select case (config_dycore)
+ case ('SIA')
+ call land_ice_sia_block_init(block, err)
case ('L1L2')
- call land_ice_lifev_block_init(mesh, err)
- case ('SIA')
- call land_ice_sia_block_init(mesh, err)
+ call land_ice_lifev_block_init(block, err)
+ case ('FO')
+ call land_ice_lifev_block_init(block, err)
+ case ('Stokes')
+ call land_ice_lifev_block_init(block, err)
case default
write(*,*) trim(config_dycore), ' is not a valid land ice dycore option.'
err = 1
@@ -278,10 +291,14 @@
err = 0
select case (config_dycore)
+ case ('SIA')
+ call land_ice_sia_solve(mesh, state, err)
case ('L1L2')
call land_ice_lifev_solve(mesh, state, err)
- case ('SIA')
- call land_ice_sia_solve(mesh, state, err)
+ case ('FO')
+ call land_ice_lifev_solve(mesh, state, err)
+ case ('Stokes')
+ call land_ice_lifev_solve(mesh, state, err)
case default
write(*,*) trim(config_dycore), ' is not a valid land ice dycore option.'
err = 1
</font>
</pre>